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SUMMARY 


A numerical method based on the conservation form of the full potential equation has 
been applied to the problem of three-dimensiofial supersonic flows with embedded sub- 
sonic regions. The governing equation is cast in a nonorthogonal coordinate system, and 
the theory of characteristics is used to accurately monitor the. type-dependent flow field. 
A conservative switching scheme is employed to transition from the supersonic marching 
procedure to a subsonic relaxation algorithm and vice versa. The newly developed com- 
puter program can handle arbitrary geometries with fuselage, canard, wing - , tbrcxvh 
nacelle, vertical tail and wake components at combined angles of attack and sideslip. Re- 
sults are obtained for a variety of configurations that include a Langley advanced fighter 
concept with fuselage centerline nacelle, Rockwell’s Advanced Tactical Fighter (ATF) with 
wing mounted nacelles, and the Shuttle Orbiter configuration. 
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1. INTRODUCTION 


An examination of the literature for supersonic/hypersonic aircraft provides an indi- 
cation of the flexibility and generality required for a prediction technique. Typical config- 
uration development variables include wing section, incidence, height, dihedral, planform, 
effectiveness of longitudinal control surfaces for trim, effectiveness of empennage for direc- 
tional stability, and propulsion system-airframe interactions. 

State-of-the-art response to these prediction requirements is provided by hypersonic 
impact methods as well as linearized analysis and design algorithms. These approaches can 
treat complex geometries with minimum response time and cost, with efficient predicted 
data coverage in terms of Mach number, angle of attack, trim deflection, yaw angle, etc. 
Shortcomings are present, however, in both the impact and linearized methods. For the 
former, interference between surface elements is totally ignored in implementations such 
as classical Newtonian, tangent wedge, and cone theories. Crossflow interactions and 
stagnation point singularities are also implicitly disregarded. In the latter, shocks, vorticity, 
and entropy wakes and layers are excluded. Furthermore, superposition of elementary 
solutions such as those for thickness and angle of attack freely used in linear models are, 
strictly speaking, invalid at hypersonic speeds. 

A need exists for new aerodynamic prediction techniques to optimize vehicles designed 
to travel at supersonic/hypersonic speeds. One requirement of a new aerodynamic pre- 
diction technique is that it be more accurate than simple noninterfering panel methods. 
Another specification is that it be more computationally efficient than currently available 
explicit finite-difference methods so that it can be incorporated into a practical design 
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procedure. The new approach should include enough of the physics of the flow to allow 
realistic optimization and should permit consideration of appropriate interactions between 
components of promising arrangements, since this has been found to be the key to in- 
creasing aerodynamic efficiency using linear methodology. Nonlinear potential theoretical 
formulations hold the promise of meeting this objective and providing economic design 
codes which are responsive to conceptual vehicle definition efforts. 

A nonlinear aerodynamic prediction technique based on the full potential equation in 
conservation form has been developed for the treatment of supersonic flows. References 1 
and 2 contain the details of the theoretical development of the full potential analysis 
method. Detailed description of the method can be found in the publications in Ap- 
pendix A. The second and third papers of the Appendix constitute a user’s manual for the 
full potential analysis code. Contained therein is a brief description of the code organiza- 
tion, main program and subroutines, flow variables, input data format, and sample test 
cases. 

The computer program entitled, “Nonlinear Supersonic Full Potential Analysis Pro- 
gram” can be obtained for a fee from: 


Computer Software Management and Information Center (COSMIC) 

112 Barrow Hall 

University of Georgia 

Athens, GA 30602 

(404) 542-3265 
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Request the program package by the designation LAR-13413. The program is written 
in FORTRAN IV for use on the Control Data 6600 and the CYBER series of computers. 
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2. FULL POTENTIAL METHOD 


The full potential equation in conservative or nonconservative form is frequently 
used for transonic flow analysis, where the local Mach number does not exceed approx- 
imately 1.4. If the assumptions of irrotationality and isentropicity are reasonably valid, 
then the full potential equation is expected to yield results comparable to Euler equations, 
even for supersonic/hypersonic flow fields. For conceptual design studies, where short 
response time is desired, the full potential methods can be an attractive substitute for 
expensive Euler methods and less accurate linear theory methods. 

A nonlinear aerodynamic prediction technique based on the full potential equation in 
conservation form has been developed for the treatment of supersonic flows. A detailed 
description of the method has been presented in several published papers 1-6 . The most 
recent publications are enclosed in Appendix A for convenience. The first three papers 6-8 
describe the method for the treatment of predominantly supersonic flows with regions of 
subsonic flow that usually occur at low supersonic Mach numberswith the 2 nd and 3 rd 
papers constituting a user’s manual for the full potential analysis code. The final two 
papers 9,10 in the Appendix describe the method to treat the unsteady form of the full 
potential equation. Additional information on the unsteady treatment may be found in 
References 11 and 12. For blunt nosed configurations with a detached bow shock, the 
unsteady method is used to generate the starting solution for the marching code. Since 
the Appendix describes the full potential method, only the results not included in the 
published articles are presented here. 
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The following summarizes the salient features of the subject full potential code. Details 
can be found in the noted references. 

• Equation in conservation form (see Refs. 1-6) 

• Flux linearized upwind differencing in the marching direction (see Refs. 2-6) 

• Conservative switch operators to treat embedded subsonic zones (see Refs. 2 and 3) 

• Treatment of wakes (see Refs. 2-6) 

• Yaw and angle of attack (see Refs. 2, 5, and 6) 

• Complex geometry treatment (fuselage, canopy, wing, canard, nacelle, tail, multibody, 
etc.) (see Refs. 2, 5, and 6) 

• Treatment of blunt nose using unsteady full potential methods (see Refs. 9-12) 

• Numerical grid generation with constraints (see Refs. 2, 4, 5, and 6-8) 

• Use of GEMPAK 13 or CDS 14 to generate geometry input files (fully automated) (see 
Refs. 2-8) 
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3. RESULTS 


The full potential analysis code can handle complex aircraft geometries as well as 
multibody configurations. The following set of four different configuration studies clearly 
demonstrate the versatility and robustness of the code in handling a wide variety of non- 
linear flows. 

The results to be presented here are 

1) Langley’s canard-wing fighter configuration, Figures 1 through 21. 

Figures 1-2 indicate the complexity of the fighter geometry. Figure 3 schematically 
shows the variation of the cross plane geometry from the nose to the back end. Figures 4- 
11 show the cross-sectional geometry and the corresponding grid setup at various axial 
stations. Figure 12 shows the circumferential pressure distribution at an axial station 
where only the canard is present. Figure 13 shows the pressure contours, cross flow velocity 
vectors, and the cross flow streamlines at this axial station. The cross flow velocity vectors 
and cross flow streamlines are obtained by projecting the total velocity vector on a unit 
sphere whose center is at the nose of the geometry. Figures 14—19 show similar results at 
other axial stations. The formation of a shock around the inlet is clearly seen in Fig. 19. 
Figure 20 shows the computational geometry and the surface grid points for this Langley 
fighter. Figure 21 shows the pressure contours in the plane of symmetry. 

2) Rockwell’s Advanced Tactical Fighter with a flow through nacelle, Figures 22 
through 24. 
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Figure 22 shows the geometry of an advanced tactical fighter with a nacelle mounted 
on the undersurface of the wing. The figure also shows the surface grid setup. Figure 23 
shows the grid at an axial station where the nacelle is present. Figure 24 shows the pressure 
contours at this axial station and the corresponding crossflow velocity 
vectors. 

3) Isolated Shuttle Orbiter flow, Figures 25 through 27 , 

Figure 25 shows the geometry of the Shuttle Orbiter. Figure 26 shows the upper 
surface chordwise pressure distribution at various span stations at = 1.4 and a = 
— 1.94°. The agreement with experimental data is very good. Figure 27 shows the OMS 
pod shock formation and its impingement on the upper wing surface. 

4) Shuttle Orbiter — External Tank mated configuration, Figures 28 through 31 . 

Figure 28 shows the multibody problem of the Shuttle mated configuration with the 
External Tank and the Solid Rocket Boosters present. Figure 29 shows a typical gridding 
for this multibody problem at an axial station where the OMS pod, External Tank/SRB, 
and the blockage are all present. Figure 2-9 shows the pressure contours at this axial station 
indicating the OMS pod shock and the detached shock in front of the blockage. Figure 30 
shows the Orbiter lower surface chordwise pressure with and without the blockage effect. 
The comparison with experimental data is good when the blockage effect is accounted 
for. Figure 3i shows the computational geometry and surface gridding for the mated 
configuration. More results on this multibody problem can be found in Ref. 15. 
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The Shuttle Orbiter and the Orbiter/External Tank studies were performed as part of 
another contract from the Space Division of Rockwell International and funded by NASA- 
MSFC, Contract No. NAS9-14000. They are included in this contract report only for 
completeness in illustrating the overall capability of the full potential code. 
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4. CONCLUSIONS 


The full potential analysis code has developed into a powerful nonlinear tool for the 
analysis of complex aerodynamic configurations: Modifications and enhancements allow 
analysis of complete configurations including fuselage, canard, wing, vertical and/or hori- 
zontal tail, nacelle (body or wing mounted), wake interference effects and multibody flows. 
Comparisons with available experimental data are in good agreement. 

The code is operational on several computer systems, such as the CRAY-XMP, CY- 
BER 203, CDC 7600, CYBER 176 and 175. A vectorized version of the code for the 
VPS-32 (modified CYBER 205) is currently in development. When fully developed and 
optimized, the vectorized code is expected to be able to analyze a complete fighter-like 
configuration in 20-30 seconds using the VPS-32 or CRAY-XMP class machine. 
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Fig. 1. Canard-wing fighter configuration. 
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Fig. 2. Treatment of complex aerodynamic configurations. 



Fig. 3. Variation in cross-sectional shape. 


15 




16 


Fig. 4. Cross-section and grid at x = 5 (fuselage/canopy/beginning of canard). 
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Fig. 5. Cross-section and grid at x = 7 (fuselage i 
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Fig. 12. Solution for Langley fighter configuration. 
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Fig. 17 . Solution for Langley fighter configuration. 
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Fig. 18 . Solution for Langley fighter configuration. 
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Fig. 22. Geometry of an advanced fighter with a nacelle. 
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Fig. 23. Grid setup for a wing mounted nacelle. 
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Fig. 24. Tactical fighter results at M 
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Fig. 30 . Orbiter lower surface chordwise pressure distribution. 
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NONLINEAR COMPUTATION OF WING-BODY-VERTICAL TAIL-WARE FLOWS AT LOW SUPERSONIC SPEEDS 


Kuo-Yen Szema* and Vijaya Shankar** 
RockweLl International Science Center 
Thousand Oaks, California 


Abstract 

A numerical method based on the conservation 
form of the full potential equation has been ap- 
plied to the problem of three-dimensional super- 
sonic flows with embedded subsonic regions. The 
governing equation is cast in a nonorthogonal 
coordinate system, and the theory of character- 
istics Is used to accurately monitor the type- 
dependent flow field. A conservative switching 
scheme Is employed to transition from the super- 
sonic marching procedure to a subsonic relaxation 
algorithm and vice versa. The newly developed 
computer program can handle arbitrary geometries 
with fuselage, wing, vertical tail and wake com- 
ponents at combined angle of attack and side- 
slip. Results are presented for a low supersonic 
mach numbers flow over the Shuttle orbiter (in- 
cluding the OMS pods and vertical tail), and for 
flows over a realistic fighter type configura- 
tion. Comparisons with experimental data are 
shown to be in good agreement for various cases. 

I. Introduction 


applicability to complex geometries and intricate 
shocked flow fields. 

The main purpose of this study is to extend 
the methodology of Reference 4 to investigate 
supersonic flows with large embedded subsonic re- 
gions over complex configurations, and as well as 
extend the treatment to combined angle of attack 
and yaw cases. All the calculations reported in 
this paper were performed using the CDC CYBER 176 
computer. A typical calculation over a complete 
configuration requires 15 minutes of CPU time on 
the CYBER system or 3 minutes on the CRAY l. 

II. Basic Formulation 

The steady, conservative full potential equa- 
tion cast in an arbitrary coordinate system de- 
fined by C * C(x, y,z), n * n(x,y,z) and ? =* 
? (x,y,z) can be written as 

r °7L + rp j \ +( ° ? 


Currently available numerical algorithms to 
compute the low supcrsoaic, inviscid flow about 
complex configurations are frequently either 
inadequate 1 or too costly to use for routine 
analysis. 2 * 3 For treatment of low supersonic 
flows, the full potential method* 4 is an ideal 
substitute for the Euler methods to avoid the 
requirements of excessive computer time and 
memory. 

Recently, Shakar et al 4 have developed a 
numerical method based on the characteristic 
theory to solve the problem of supersonic flow 
with embedded subsonic regions. Reference 4 
describes the characteristic theory involved in 
determining the condition for a marching direction 
to exist. Once that condition is violated, the 
marching scheme Is transitioned to a relaxation 
scheme through a conservative switching opera- 
tor. For the marching condition violation, the 
total velocity q does not have to be subsonic. 

Even for a supersonic total velocity q, if the 
component in the marching direction is subsonic, a 
relaxation scheme is required. In order to pro- 
perly produce the necessary artificial viscosity 
through density biasing, Reference 4 defines two 
situations: (L) the total velocity q is super- 

sonic, but the marching direction component is 
subsonic (defined as Marching Subsonic Region 
(MSR) ) , and (2) the total velocity q Is subsonic 
(termed as Total Subsonic Region (TSR)). The 
method of Reference 4 uses a numerical mapping 
technique to generate the body fitted, nonortho- 
gonal curvilinear coordinates system. The key 
advantage is that it lias no restrictions on its 


where the density p is given by 


0 - [l - ( ^Y~ 1 M* (u*„ + VJ^ + W> p - ll] 

(2) 


and is the free stream mach number, a is the 
local speed of sound and U,V,W are the contravari- 
ant velocity components. Introducing the follow- 
ing notation for convenience. 

Ui = U U 2 * V Uq « W 

x » XI y * X2 z = xq 

r - X l n - X2 p 3 x 3 


the cont ravariant velocity can be expressed as 
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The Jacobian of the transformation J is 
represented by 
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The nature of Eq. (1) can be analyzed by 
studying the eigenvalue system of Eq. (1) combined 
with the irrotationality condition in the (C,n) 
and (C ,F) planes. A detailed discussion on this 
can be found in Reference 4, Therefore, only the 
final results are presented here. 
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At a grid point, the marching direction K 
is hyperbolic and the total velocity q is 


supersonic, fa - — 1 <0, q > a. This 
11 a 2 

point will use the algorithm of Reference 

5. 


2. At a grid point, the marching direction C 

U 2 > 


is elliptic, f a 
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-) > 0, but the 


total velocity q is supersonic, q > a. 
This point will be treated by a transonic 
operator with a built-in density biasing 

, a 2 

based on the magnitude of (1 1. 

r\2 


3. At a grid point, the direction C is 

elliptic and the total velocity q is sub- 
sonic , q < a. This point will be treated 
by a subsonic central differenced 
operator. 


III. Method of Solution 


Figure 1 shows the schematic of a fuselage- 
canopy forebody geometry with an embedded MSR and 
TSR present in a supersonic flow. To solve this 
problem, the marching scheme of Reference 5 is 


used when f < 


U 2 


11 


] is negative, and a relaxation 


s jj2 

scheme is used when fa ni ) is positive. 

11 7 

a* 

First, march from the nose up to the plane denoted 
by (A-B) in Fig. 1, using the method of Reference 
5. Then, between (A-B) and (C-D), which embed the 
subsonic bubble (MSR and TSR) , use a relaxation 
scheme and iterate until the subsonic bubble is 
fully captured. Then, resume the marching scheme 
from the plane (C-D), downstream of the body. 
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Fig. 1 Embedded subsonic bubble in a supersonic 
flow. 


1. Treatment of ft /AC (p U/J) Term 

The finite-difference operator at point 
(i+1, ft,k) for the first terra in Eq. (1) may be 
written as 


ft 

AC 


n r U> & t U\ 

— i P T 1 * 9 i ' D “ 


J i+1 


supersonic 

• + 

+ U - 9 i+i> Kc ^ j\+i ^ 
marching 

where subsonic 

ft refers to backward differencing 
3 refers to forward differencing 

■ 1 if fa n - < 0 

a 1 

t U 2 

• 0 if (a„ - -> 0 
11 2 

The first terra in Eq. (4) corresponds to the 
supersonic marching operator and the second term 
is the subsonic operator. By using a local iiner- 
ization procedure, Eq. (5) can be expressed in 
terra of <b only. Details of the procedure are 
given in Refs. 4 and 5. 


2. Treatment of ft/ftti (p V/J) Term 


The finite difference operator for the second 
terra in Eq. (1) is given as 


ft_ r V> 

ftT) P J 


— (~ h 

1+1 ° J J+1/2 
supersonic 


< 6 > 

where marching 

^2 subsonic 

A *1 if fa - — 1 < 0 (supersonic 

a 2 i+1 point) 


6,,, - 0 if fa. - — I >0 (MSR) 
1+1 11 a 2 i+1 


When * 1» that is, the point is super- 

sonic with respect to C , only the first_term in 
Eq. (6) is used and the biased density p is 
defined by (for V > 0) 


°j+l/2 (1 v j+l/2 )p j+1/2 + 2 v j+l/2 (D j + 

where 

- a 2 

v * maxf 0,1- a — ) 

22 V 2 

In Eq. (7), the evaluation of o* depends on 
whether the flow is conical or nonconical. For 
conical flows, all o* quantities are evaluated at 
the i th plane. For nonconical flows, at each non- 
conical marching plane, initially p* is seen to be 
the value at the i c ^ plane and then subsequently 
iterated to convergence by setting p* to the 
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previous iterated value of p at the current i+l 
plane. 

When the point is e lliptic , the density 
biasing is defined by 


0 j+l/2 


where “v * raaxf 0, 



As before, the super- 


script n+1 denotes the current relaxation cycle 
for a subsonic bubble calculation. ^Note the dif- 
ference in the definition of v and v'. The density 
biasing in the cross flow direction n is turned 
off when the total velocity q is less than the 
speed of sound a. The implicit treatment of V in 
the marching subsonic operator of Eq (6) is the 
same as that of the supersonic part, explained in 
Reference 5* 


r 

A similar procedure is implemented for ip yip, 
terra in Eq. (1). 

3. Implicit Factorization Algorithm 

Combining the various terms of Eq. (1) as rep- 
resented by Eqs. (5)-(8) together with the terms 

r W> 

arising from fp will result in a fully ira- 


The density p appearing in Eq. (9) and Eq. (10) 
can be either p or p depending on the sign of 


a ll as illustrated in Eq. (6). 

a 2 

Equation (9) has the form L^L (A6) 
and it is implemented as follows: 11 


U(Ad>)* =» R L (A<b ) » (AA )* <fc * A , + Ait* 
' n i+L i 


The various quantities appearing in Eq. (9) are 
given by 

fe 
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plicit model. This is solved using an approximate 
factorization implicit scheme. After some rear- 
rangement of the terms, the factored implicit 
scheme becomes 
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_l A pa 22 A 
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n - <? ) A t_ (“Jlii'i 
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Ar o -h+2 - r i+i 


A r * r - r 
i+t i 


(ID 


If the flow field does not contain an embedded 
MSR or TSR, the implicit factored algorithm of 
Eq. (9) performs a pure marching procedure start- 
ing from an initial known data plane. In this 
situation, there is no need to go back to the 
upstream starting plane and iterate the solu- 
tion. However, if a subsonic bubble is presnt 
(between planes AB and CD in Fig. 1, then the 
solution procedure of Eq. (9) performs a relaxa- 
tion method, and iterates for the elliptic sub- 
sonic bubble to converge. 

IV. Initial Conditions 


For a pure supersonic flow, initial conditions 
are required at the starting plane. Usually, the 
starting plane Is set close to the nose region of 
the configuration. For a sharp nosed configura- 
tion, conical solutions are prescribed, and for a 
blunt nose, the axisymmetric unsteady full poten- 
tial solver of Reference 6 is used to obtain flow 
field in the transonic forbody region. 


I 2_ 

r dr. 


k v ♦ 


pa 


_22 5- * 

J AF *i 


1 <°h? 1-1 A* 1 !_ 2_\ aa 

A An J af i+l A *i+l A AF J An i+t AP i+l 

( 10 ) 


In the embedded subsonic region, when Eq. (5) 
is applied at an i+l grid point, information on 
the flux oU at i+2 is required. For the first 
relaxation pass, sonic conditions are assumed at 
i+2. 
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°i+2 


D* 



m 2 1 l/y-l 

as 


„ * / xl/2 

U i+2 = q (a l 1 J i+2 


where 


!d* y_1 /m 2 ] 1/2 


( 12 ) 


The sonic values p* and q* are purely a func- 
tion of the free-stream Mach number M^. For the 
second relaxation cycle and onwards, the condi- 
tions from the previous relaxation cycle is used. 


V. Boundary Conditions 


In order to solve the full potential equation, 
it is essential to specify appropriate boundary 
conditions of the body surface and outer boundary. 

1. Body Surface 

At a solid boundary, the contravariant veloc- 
ity V is set to zero. Exact implementation of V = 
0 in the implicit treatment of Eq. (9) is des- 
cribed in Reference 4. 

2. Outer Boundary 

The outer boundary is set away from the bow 
shock and the freestream velocity potential <t> m is 
imposed along that boundary. All discontinuities 
in the flow field are captured. The precise den- 
sity biasing activator v, based on the character- 
istic theory, allows for sharp capturing of shocks 
in the flow. 

3. Symmetric Boundary Conditions 

For yaw angle 8*0, only the half plane prob- 
lem needs to be solved with the plane of symmetry 
boundary conditions imposed along K ■ 2 and (Kmax 
- 1), as shown in Fig. 2a. Imposing that the flow 
conditions along K * 1 are to be the same as the 
ones along K * 3, the L^ operator results in a 
tridiagonal system that can be easily solved. 
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a) PLANE OF SYMMETRY 
(YAW ANGLE 0) 

"i.i ,;, i*3 

■j.KMAX KMAX-2 
Fig* - Symmetry and 


KMAX 1 



b) PERIODIC CONDITION FOR 
YAW TREATMENT 

‘‘'i.i = ^j.KMAX- 2 

<-’j, KMAX ' '-V3 

periodic conditions. 


4. Combined Yaw and Angle of Attack 


Even for a symmetric configuration, when yaw 
angle is present the entire cross-flow plane needs 
to be solved as shown in Fig. 2b. In this case 
the flow conditons along K * 1 are set to be the 
same as the ones along K * (KMAX - 2). This 
destroys the tridiagonal nature of the 
operator. A special routine has been developed to 
invert a matrix of the following type. 


XX X 0 

XXX 00 

XXX 0 


(13) 


0 0 XXX 

OX XX 


In the current formulation, positive angle of 
attack a represents a positive cartesian velocity 
v in the freestream and similarly positive yaw ft 
produces a positive w in the free stream. When 
both angle of attack and yaw are present, first 
the freestream is turned by an angle R and then by 
a. 


Let (x,y,z) be the inertial Cartesian sys- 
tem. After an initial yaw turn 8 let the wind 
axis system be (x* ,y f ,z f 2,^and after a subsequent 
a turn let it become (x,y,z). 


X 

cosa sina 0 


y = 

-sina cosa 0 


z 

o 

o 



cosa cosR 

sina 

= 

-sina cosfi 

cosa 


-sin8 

0 



The free stream is now along x. The normalized 
free stream velocity potential is given by 


X cosa cosB + y sina + z cosa sinft 


(15) 


Using Eq. (14), the lift, drag and side forces are 
easily represented. 


F cosa cosR + F sina + F cosa sinP 
x y z 


L * — F sina cosfi + F cosa (16) 

x y 

5. Swept Trailing Edge Wake Treatment 

Figure 3 shows a schematic of a swept trailing 
edge wake system. In order to treat the region 
behind the trailing edge, an artificial cut is 
created and the pressure jump [P ] across this cut 
is imposed to be zero as a boundary condition. 

This is achieved by maintaining the jump in the 
velocity potential <t> along a k * constant line 
(see Fig. 3) for j * 2 to be the same as the value 
[<t>] at the trailing edge. The full potential 
equation is not solved at grid points on the wake 
cut. Instead, * 0 is solved to provide [i» ] = 

0 across the wake cut. Maintaining [6] constant 
along a k line provides [ d>r ] * 0. The combination 
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SOLVE 0 AT WAKE POINTS 

Fig, 3 Wake boundary condition. 


of [<b r ] = 0 and [d^] * 0 across the cut satisfies 
[?] *’0 approximately. 


6. Geometry and Grid System 


• KEY POINTS ATINPUT 
GEOMETRY STATONS (x a. c) 

A INTERPOLATED KEY POINTS 
AT A MARCHING PLANE (x b) 

y V I. II. Ill PATCHES 



Fig. 4 Geometry setup. 


Figure 5 shows the pressure distribution on the 
surface of an arrow-wing configuration at location 
x/S =0,8 for M m = 2.9, and a = 10°. The improve- 
ment in the prediction capability achieved using 
the wake treatment is illustrated. The dashed 
line represents the result from "no wake" treat- 
ment (assumes a flat plate behind the trailing 
edge) and the solid line represents the modifica- 
tions to the pressure distribution once a zero 
jump in pressure across the wake cut is imposed. 
The soLid line pressures on the body agree very 
well with experiments. Without a proper wake cut 
treatment, the overall lift and drag forces and 
the pitching moment can be off by a considerable 
margin. 


The geometry of a configuration is prescribed 
at discrete points in a crossplane (usually x = 
constant plane) at various axiaL locations. These 
geometry input points are usually obtained from a 
geometry package such as GEMPACK or CDS. The in- 
put points are then divided into several patches, 
and at each patch a key-point system is estab- 
lished as shown in Fig. 4. The geometry at a 
marching plane is then obtained by joining the ap- 
propriate key-point for each patch. Using a cubic 
spline passing through the key points, a desired 
grid point distribution (clustering) is set up on 
the body surface. Then, using an appropriate 
outer boundary, the grid for the flowfield calcu- 
lation is generated by using an elliptic grid 
generator. 

VI. Results and Discussion 

Four cases are presented to substantiate the 
recently developed code. 



1. Flow over an arrow-wing body at ^ 3 
2.96, a = 10.01°. 

2. Flow over a forebody configuration at 

= 2.5, R = 5° and at = 1.7, 

a = 10°, 3 = 5°. 

3. Flow over the entire shuttle orbiter 
geometry at = 1.4, a * 0°. 

4. Flow over a realistic fighter configura- 
tion at different angles of attack and 
freestream mach numbers. 


Fig. 5 Circumferential pressure distribution 
for the arrow-wing at x/9 = 0.8L, 

= 2.96, a = 10°. 

Figure 6 presents the pressure distribution on 
a fully developed forebody® for >^=2.5 and yaw 
angle 3 = 5°. Figure 7 shows the circumferential 
pressure distribution for the same body at k/9 = 
0.68 for = 1.7, or = 10° and Q » 5°. The exper- 
imental data 7 are also given in these two figures. 
The results show that the present prediction is in 
excellent agreement with the experimental data. 
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Fig. 6 Pressure distribution on a Eorebody in 
sideslip * 2.5, 3 - 5.02°, a - 0°. 
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Fig . 7 Circumferential pressure distribution on 
a forebody at x/£ = 5.02°, a ■ 0°. 

Figures 8 to 18 give the geometry and the cor- 
responding flow field solutions for the shuttle 
orbiter. The side view, cross section and grid in 
the nose region are given in Fig. 8. The top 
views, cross-section, grid and chordwise cross- 
section are presented in Fig. 9. It should be 
mentioned here that for most of the previous space 


SC83 25232 



Fig. 9 Geometry and grid generation for Shuttle 
orbiter at axial cuts. 

shuttle calculations 9 » 10 the geometry has been 
modified by smoothing out the canopy and increas- 
ing the wing sweep angle from 45 to 55 degrees in 
order to avoid the subsonic bubble. Since the 
present method Is valid for supersonic flow with 
embedded subsonic region, the realistic shuttle 
orbiter geometry was used without any 
modification. 

Figure 10 shows the surface pressure distribu- 
tion along the leeward plane of symmetry. At x 5 
170 in. which is the beginning of the canopy, the 
pressure increases rapidly from Cp = 0.3 to 1.0, 
approximately. In the canopy region an MSR/TSR is 
formed and required three relaxation cycles to 
develop the solution. The circumferential surface 
pressure distribution at x » 240 in. is shown in 
Fig. 11. The experimental data is also given in 
these figures and the agreement is very good. The 
surface pressure distribution along the wing lead- 
ing edge is given in Fig. 12. It is seen that the 
present calculation agrees with the experiment 
data quite well along the entire wing leading edge 
except in the vicinity of the wing-fuselage junc- 
tion, where a vortex flow may exist. 

Figure 13 to Figure 17 present the orbiter 
chordwise pressure distribution on the upper and 
lower wing surface at z * z/bw = 0.471, 0.530, 
0.641, 0.780 and 0.887, respectively where bw is a 
semispan defined in Fig. 8. The results show that 



51 



SCI32I314 



Fig. 10 Surface pressure distribution at leeward 
plane of symmetry- 


Ktl 1*1* 

1 .5, 1 F T ! ! 1 I I | 



-0.51 1 I 1 1 i I I 1 * 

0 20 40 60 80 100 120 140 160 180 

ft 

Fig* 11 Surface pressure distribution around the 
body. 
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Fig. 12 Surface pressure distribution along the 
wing leading edge. 

the present predictions are in very good agreement 
with the experimental data, except in the region 
near the trailing edge of the upper surface at 
z - 0.471 and 0.53 span stations. Here again, a 
vortex flow or separation may be causing the dis- 
crepancies. Figure 18 shows the circumferential 
pressure distribution for the orbiter at x = 1120 
in. It is noted that the pressure at the vertical 
tail and OMS pods are well predicted. 


ten JH»» 



Fig. 13 Shuttle orbiter chordwise pressure 
distribution; 3 1.4, a 3 0.0°, 
z * 0.471. 



Fig. 14 Shuttle orbiter chordwise pressure 
distribution; ^ ■ 1.4, a * 0.0°, 
z » 0. 53. 
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Fig. 15 Shuttle orbiter chordwise pressure 

distribution; M*> = ^ ^ a * 0.0°, 
z = 0.641. 
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Fig* 16 Shuttle orbiter chordwise pressure 
distribution; * 1.4, a * 0.0°, 
z - 0.78. 
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Fig. 17 Shuttle orbiter chordwise pressure 
distribution; =■ 1.4, a = 0.0°, 
z = 0.887. 


Figure 20 presents the surface pressure at 
various axial stations and the corresponding grid 
distribution for the wing body geometry. For this 
case, the MSR/TSR starts around x - 0.4 at the 
leading edge, and remains subsonic all the way to 
the end of the wing. Figure 21 and 22 show the 
circumferential pressure distribution in the 
vertical tail and wake region of the fighter-like 
configuration. The pressure on the vertical tail 
surface is given separately along y-direction. 

The results clearly show that the present wake 
tretraent provides the correct zero pressure jump 
condition across the wake. The chordwise pressure 
distributions from the center of the body to the 
top of the wing are given in Fig. 23. Figure 24 
presents the pressure contour on the upper and 
lower surface of the fighter configuration. 

The lift and drag coefficients from the 
present calculations for this fighter model are 
also given in Table l. The comparison with 
experimental data show excellent agreement. 

VII. Conclusions 

A nonlinear full potential method has been 
applied to investigate the supersonic flows with 
embedded sobsonic regions over some very complex 
configurations. A* conservative switching scheme 
is employed to transition from the supersonic 
marching algorithm to a subsonic relaxation pro- 
cedure. The present predictions are in very good 
agreement with experiment data. 

Work is now progressing to simulate the multi- 
body interaction between the shuttle orbiter and 
the external tank. The present methodology will 
also be extended to treat all mach number flows 
(fully subsonic flow as well as subsonic flow with 
pockets of supersonic region (transonic case)). 
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I. THEORY 


ABSTRACT 

A numerical method based on the conservation form of the full poten- 
tial equation has been applied to three-dimensional supersonic flows with em- 
bedded subsonic regions. The governing equation is cast in a nonorthogonal 
coordinate system, and the theory of characteristics is used to accurately 
monitor the type -dependent flow field. A conservative switching scheme is 
employed to transition from the supersonic marching procedure to a subsonic 
relaxation algorithm and vice versa. The newly developed computer program can 
handle arbitrary geometries with fuselage, wing, vertical tail and wake com- 
ponents at combined angles of attack and sideslip. Example results in this 
report include the Shuttle Orbiter flow at a low supersonic Mach number, flow 
over a realistic fighter-type configuration, wake simulations for an arrow 
wing, and a forebody in sideslip. Comparisons with experimental data are 
shown to be in good agreement for various cases. 

INTRODUCTION 

Nonlinear aerodynamic prediction methods based on the full potential 

\ 2 3 5 

equation are used regularly for treating transonic ’ and supersonic flows 
over realistic wing-body configurations. The transonic algorithms 1 are de- 
signed to treat predominantly subsonic flows with pockets of supersonic re- 
gions bounded by sonic lines and shocks. The supersonic methods 3-5 are based 
on a marching concept, and require the flow to remain supersonic in a given 
marching di recti on* Once the marching direction velocity becomes subsonic, 
the domain of dependence changes and a pure marching scheme will violate 
the rules of characteristics signal propagation. The possibility of a march- 
ing velocity becoming subsonic in a supersonic flow is great, especially for 
low supersonic freestream Mach number flows (M w = 1.3 - 1.7) over moderately 
swept fighter-like configurations (sweep angle A 45 - 50°) and over forebody 
shapes having a sizeable fuselage-canopy junction region. The methodology of 
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Refs. 6 and 7 is an extension to the marching scheme of Ref. 5 and is designed 
to treat embedded subsonic regions in a supersonic flow. In order to properly 
produce the necessary artificial viscosity through density biasing, Ref. 6 de- 
fines two situations: 1) the total velocity, q, is supersonic, but the march- 

ing direction component is subsonic (defined as Marching Subsonic Region 
(MRS)), and 2) the total velocity, q, is subaonic (termed as Total Subsonic 
Region (TSR) . 

The method of Refs. 4-7 is based on the characteri Stic theory of sig- 
nal propagation and uses a generalized, nonorthogonal , curvilinear coordinate 
system. It has no restrictions (limitations of the full potential theory 
hold) on its applicability to complex geometries and intricate shocked flow 
fields. It is a conservative formulation and uses numerical mapping tech- 
niques to generate the body-fitted system. 

This report presents a brief description of the overall methodol - 
4-7 

ogy along with some user information on code organization, input and output 
data, and sample results. A typical calculation over a complete configuration 
(wing-body-vertical tail -wake) requires fifteen minutes of CPU time on the 
Cyber 176 machine or three minutes on the Cray-1. 

BASIC FORMULATION 

The steady, conservative full potential equation cast in an arbitrary 
coordinate system defined by c * C(x,y,z), n = n(x,y,z) and r = r.(x,y,z) can 
be written as 

( p j'e + ( p i\ + ( p tV 0 (i) 

where the density p is given by 

P * Cl - Mf fU* c + V0> n + W* K - 1}] (2) 

and is the free stream Mach number, a is the local speed of sound and U,V,W 
are the contravariant velocity components. Introducing the following notation 
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for convenience 


Ul 5 U U 2 » V 
x * xi y » x 2 
C * Xi Vi • X 2 


U3 = w 
2 * X 3 

l * x 3 


the contravariant velocity can 


be expressed as 


3 

U, * l a..d> i * 1,2,3 

j*l 1J X j 



3 ax. dXj 

k*l dx k ax k 


i » 1,2,3 (transforma- 
j * 1,2,3 tion metrics) 


(3) 


The Jacobian of the transformation J is represented by 


*(x,y,z) 


T. 

Z 


c 


z 


n 


z 


l 


1 


(4) 


The nature of Eq. (1) can be analyzed by studying the eigenvalue sys 
tern of Eq. (1) combined with the i rrotati onal i ty condition in the (c,n) and 
(C,F,) planes. A detailed discussion on this can be found in Reference 6. 
Therefore, only the final results are presented here. 


1 . 


At a grid point, the marching direction C is hyperbol ic and the total 
velocity q is supersonic, (a., - —1 < 0, q > a. This point will use 
the algorithm of Reference 5. 
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2 


At a grid point, the marching direction c is elliptic, fa,, - —'I > 0, 

II 2 

but the total velocity q is supersonic, q > a. This point wilt be 
treated by a transonic operator with a built-in density biasing based 
on the magnitude of (1 - —1. 

q 2 

3. At a grid point, the direction c is elliptic and the total velocity q 

is subsonic , q < a. This point will be treated by a subsonic central 
differenced operator. 


METHOD OF SOLUTION 


Figure 1 shows the schematic of a fuselage-canopy forebody geometry 
with an embedded MSR and TSR present in a supersonic flow. To solve this 
problem, the marching scheme of Reference 5 is used when fa,, - — 1 is nega- 

i|2 ' l 2 

tive, and a relaxation scheme is used when fa,, - — 1 is positivi. First, 
march from the nose up to the plane denoted by (A-§f in Fig. 1, using the 
method of Reference 5. Then, between (A-B) and (C-D) , which embed the sub- 
sonic bubble (MSR and TSR), use a relaxation scheme and iterate until the sub- 
sonic bubble is fully captured. Then, resume the marching scheme from the 
plane (C-0), downstream of the body. 


Treatment of a/ac (o U/J) Term 

At a grid point (i + 1, j , k ) , the derivative in the marching direc- 
tion c is given by 


o ( u, 
dC T 


9 i k fu 1 


where 


supersonic 


+ (1 - ®i + i) ac (o 

marching 

subsonic 


( 5 ) 


a refers to backward differencing 


6l 


5 refers to forward differencing 


04 ■ 1 If (a,, - — ' ) < 0 
i 11 a 2 

y2 

» o if (a,, - — > 0 
11 a 2 

The first term in Eq. (5) corresponds to the supersonic marching 
operator and the second term is the subsonic operator. By using a local lin- 
earization procedure, Eq. (5) can be expressed in terms of * only. Details of 
the procedure are given in Refs. 4 and 5. 


Treatment of S/S-n (p V/J) Term 


tn fo J 1 ' 6 i+l 5n ,0 iVl/2 




supersonic 


= V> 




1/2 


( 6 ) 


marching 

subsonic 


When e i+1 = 1, that is, the point is supersonic with respect to _, 
only the first term in Eq. (6) is used and the biased density o’ is defined by 
(for V > 0) 


p j+l/2 = (1 " v j+l/2 )p j+l/2 + 2 v j + l/2 (p j + p j-l^ 

a 2 N 


where v = max(0, 1 - a 


22 


V 2 


( 7 ) 


In Eq. (7), the evaluation of p* depends on whether the flow is coni- 
cal or nonconi cal. For conical flows, all p* quantities are evaluated at the 
plane. For nonconical flows, at each nonconical marching plane, initially 
p* is set to be the value at the 1*" plane and then subsequently iterated to 
convergence by setting p* to the previous iterated value of p at the current 
i+1 plane. 
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When the point is elliptic in the marching direction, the density 
biasing is defined by 

*n+l M * x n , 1 » / n . n , 

p j+l/2 * ^ " v J+l/2^°j+l/2 2 \j+l/2 ^ p j p j-l^ 

2 

where v * max(0, 1 As before, the superscript n+1 denotes the current 

relaxation cycle for i subsonic bubble calculation. Note the difference in 
the definition of 7 and v. The density biasing in the cross-flow direction n 
is turned off when the total velocity, q, is less than the speed of sound a. 
The implicit treatment of V in the marching subsonic operator of Eq. (6) is 
the same as that of the supersonic part, explained in Reference 5. 

w 

A similar procedure is implemented for fp -jl^ term in Eq. (1). 
Implicit Factorization Algorithm 

Combining the various terms of Eq. (1) as represented by Eqs. (5)-(8) 

w 

together with the terms arising from fp will result in a fully implicit 
model. This is solved using an approximate factori zation implicit scheme. 
After some rearrangement of the terms, the factored implicit scheme becomes 



1 ft pa 22 ft , A . 

? T” ft7 ] 


( 9 ) 


The density p appearing in Eq. (9) can be either p or p depending on the sign 
of fa^ -—1 as illustrated in Eq. (6). 


fol 1 ows : 


Equation (9) has the form L^L (a*) « R and it is implemented as 
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L-(A4>)* = R, L (A<t>) * <b. , * 4>, + A6 

z* T| I +1 1 


INITIAL CONDITIONS 


For a pure supersonic flow, initial conditions are required at the 
starting plane. Usually, the starting plane is set close to the nose region 
of the configuration. For a sharp nosed configuration, conical solutions are 
prescribed, and for a blunt nose, the axi symmetric unsteady full potential 
solver of Reference 8 is used to obtain flow field in the transonic forebody 
region. 

In the embedded subsonic region, when Eq. (5) is applied at an i + 1 
grid point, information on the flux pU at i+2 is required. For the first re- 
laxation pass, sonic conditions are assumed at i+2. 


i+2 


•_L_ + nl 

•Y+l Y+l 


M 2 ! 1 ^ 

00 


u 


i+2 


a* (a ) l/Z 
q la llM+2 


( 10 ) 


where 


r *Y~1 /u 2-,l/2 
q* * [p* r /H m j 

The sonic values p* and q* are purely a function of the free-stream 
Mach number M«. For the second relaxation cycle and onwards, the conditions 
from the previous relaxation cycle is used. 

BOUNDARY CONDITIONS 

In order to solve the full potential equation, it is essential to 
specify appropriate boundary conditions on the body surface and at the outer 
boundary. 
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1. 


Body Surface 


At a solid boundary, the contravariant velocity V is set to zero. 
Exact implementation of V * 0 in the implicit treatment of Eq. (9) is des- 

fr 

cribed in Reference 4. 

2. Outer Boundary 

The outer boundary is set away from the bow shock and the free -stream 
velocity potential is imposed along that boundary. All discontinuities in 
the flow field are captured. The precise density biasing activator, v, based 
on the characteri Stic theory, allows for sharp capturing of shocks in the 
fl ow. 


3. Symmetric Boundary Conditions 

For yaw angle 0*0, only the half plane problem needs to be solved 
with the plane of symmetry boundary conditions imposed along K = 2 and (KMAX - 
1), as shown in Fig. 2a. Imposing that the flow conditions along K = 1 are to 
be the same as the ones along K = 3, the l P operator results in a tri diagonal 
system that can be easily solved. 


4. Combined Yaw and Angle of Attack 

Even for a symmetric configuration, when yaw angle is present the en- 
tire cross-flow plane needs to be solved as shown in Fig. 2b. In this case 
the flow conditions along K * 1 are set to be the same as the ones along K = 
(KMAX - 2). This destroys the tridiagonal nature of the l r operator. A spe- 
cial routine has been developed to invert a matrix of the following type. 


L 


T. 


XX X 0 

XXX 00 
XXX 0 


0 0 XXX 

O X XX 


( 11 ) 
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In the current formulation, positive angle of attack a represents a 
positive Cartesian velocity v in the free-stream and similarly positive yaw a 
produces a positive w in the free-stream. When both angle of attack and yaw 
are present, first the free-stream is turned by an angle 0 and then by a. 

The normalized free stream velocity potential is given by 

4 > * x cosa cose + y sina + z cosa si na (12) 

© 

The lift and drag forces are represented by 

0 * cosa C0S8 + Fy sina + ? z cosa sine 

L * -F sina cose + F cos« (13) 

x y 

5. Swept Trailing Edge Wake Treatment 

Figure 3 shows a schematic of a swept trailing edge wake system. In 
order to treat the region behind the trailing edge, an artificial cut is cre- 
ated and the pressure jump [p] across this cut is imposed to be zero as a 
boundary condition. This is achieved by maintaining the jump in the velocity 
potential <t> along a' k » constant line (see Fig. 3) for j * 2 to be the same as 
the value [<t>] at the trailing edge. The full potential equation is not solved 
at grid points on the wake cut. Instead, o » 0 is solved to provide [o n ] * 

0 across the wake cut. Maintaining [<s] constant along a k line provides 
[$£ ] * 0. The combination of [<j> c ] * 0 and * 0 across the cut satisfies 
[pj » 0 approximately. 

6. Geometry and Grid System 

The geometry of a configuration is prescribed at discrete points in a 
crossplane (usually x * constant plane) at various axial locations. These 
geometry input points are usually obtained from a geometry package such as the 
GEMPACK or COS. The input points are then divided into several patches, and 
at each patch a key-point system is established as shown in Fig. 4. The geom- 
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etry at a marching plane is then obtained by joining the appropriate key-point 
for each patch, as shown in Fig. 5. Using a cubic spline passing through the 
key points, a desired grid point distribution (clustering) is set up on the 
body surface. Then, using an appropriate .outer boundary, the grid for the 
flow-field calculation is generated by using an elliptic grid generator. ^ 

II. CODE ORGANIZATION 

The program is written in FORTRAN language. It can be executed on 
any CDC machine (Cyber 176, CDC 7600), as well as on the Cray 1. At present, 
the code is not optimized for a vector machine like the Cray or Cyber 205. 
for a cross-plane (n*C) grid of 30 x 60, the program requires a storage of 
230,000 words octal. The program consists of a main routine and several sub- 
routines. A brief description of the main program and other pertinent subrou- 
tines are given in this section. 

Program Main 

Program Main coordinates the entire operation. A flowchart describ- 
ing the various operations performed by the Main program is given in Fig. 6. 
The Main program sets up the initial (known) data plane and the body-fitted 
grid system, and performs the and operators to advance the solution. 

The marching step size AC can either be prescribed or computed at each march- 
ing plane from a given Courant number and the maximum eigenvalue. The various 
read and write tapes used in the calculation are listed below. 

Program Main (Tape 1, Tape 2, Tape 3, Tape 4, Tape 5, Tape 7, 

Tape 8, Output, Tape 6 * Output) 

Tape 1: Output solutions for plot. 

Tape 2: Output solutions for restart. 

Tape 4: Output solutions for restart. 

Tape 3: Read in starting solutions. 

Tape 5: Input data. 
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Tape 6: Solution output. 

Tape 7: Read tape containing solutions for subsonic region. 

Tape 8: Write tape for subsonic bubble calculation. 

Subroutine EIGEN (EIGENY, EIGENZ) 

This subroutine computes the maximum eigenvalue EIEGNY in the (C.n) 
plane and the maximum eigen value EIGENZ in the plane. The expression 

used for calculati ng^the eigenvalue is given in Ref. 5. The maximum eigen- 
value information is then used to compute a marching step-size AC for a speci- 
fied Courant number. 

Subroutine NFORCE (PX, PY, PM, AREA. KG) 

At the end of each marching plane calculation, this subroutine com- 
putes the axial force, PX, vertical force, PY, and the side force, PM, by 
integrating the pressure force acting on an elemental area, dA. The elemental 
area, dA, is computed from the transformation matrice using the formula (at a 
body point j ■ 2) 

dA = {[y^-z^] 2 + [x ^z ; -x ^ ^ + Cy^-y^] 2 } 1/2 dCdC . 

KG * 0, conical or blunt body nose force calculation 
* 1, rest of the body force calculation. 

The program also prints the C L and Cq information based on a prescribed refer- 
ence area, and C M about a given reference point (X 0 , Y 0 ). 

Subroutine GEOM (N9. NRP) 

N9 * 0, geometry data at X^ are read in 

* 1, geometry data at X 1 are updated 

NRP = 0, X-plane geometry calculation 

* 1, spherical plane geometry calculation. 
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Subroutine GEOM sets up the body grid points from a prescribed geom- 
etry shape. As illustrated in Fig. 4, the geometry is input in different 
patches (the number of patches to be used is left to the discretion of the 
user). From the input geometry points, a key point system is established 
using cubic splines. These key points are then joined from one prescribed 
geometry station to the other to provide the geometry at any intermediate 
marching plane. A flowchart describing GEOM is given in Fig. 7. 

Subroutine Grid 

Once the body points are obtained at a marching plane from GEOM, sub 
routine GRID sets- up the entire crossflow plane grid using an elliptic grid 
solver that satisfies certain grid constraints. Figure 8 gives the flowchart 
for GRID. 


Subroutine Metric 

This subroutine computes all the necessary transformati on metrics and 
Jacobi ans at various node and half node locations as required by the solution 
algorithm (L s and l operators). 

Subroutine UVW 

This subroutine computes all the contravariant velocities U, V and W, 
and the density p. 

Subroutine RHOBIAS 

This subroutine performs the density biasing in the (n.S) plane based 
on a characteristic theory. This operation is essential to treat crossflow 
supersonic regions and to capture shock waves. 

III. INPUT DATA 

Input data includes specifications of flow parameters, grid param- 
eters, read and write tape parameters, and geometry data in patches at various 
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axial stations. A sample inut is presented here. 

Input cards 1 through 45 are sel f-ex plained. The rest of the input 
cards are explained below. 



Format 

Symbol 

Descri ption 

Card 46 

Col. 1-80 

10F8.4 

ZTAPT(IO) 

Location of detail flow field 
printout 

Card 47 

Col . 1-5 

15 

ISC 

Number of patch (geometry) 

Card 48 

Col . 1-50 

1015 

NPT (10) 

Number of output points in 
each patch 

Card 49 

Col. 1-15 

F15.6 

XI 

x location at which th'e geometry 
cross-section is given. 

Col. 16-20 

15 

ISC1 

Number of input patch at this 
location 

Card 50 

Col . 1-5 

15 

ITH 

Patch number 

Col. 6-10 

15 

IPT 

Number of input points 

Col. 11-15 

15 

NO 

Clustering control parameter 
0: equal space 

1: clustering at the beginning 
of the patch 

2: clustering at the end of patch 

Card 51 

Col . l-i5 

F15.6 

y 

y-val ue 

Col. 16-30 

F15.6 

z 

z-val ue 


70 


120 NMAX 15 NO. OF STREAMUISE STEP 












• 

z 

























*- 

o 

























3 

44 

























X 

H- 

























1- 

<t 

























3 

X 









• 
















O 

Ui 









Ui 









cu 








1- 








• 

X 









/N 





• 


A 

44 

* 







Ui 

X 









c 





in 


pJ 









1- 

K 









(X 


M 



Z 


Ui 

A 




A 




« 

<r 



• 






z 





o 


44 

o 




c 




X 

X 




• 



• 




> 



►4 


u» 

•J 




<c 




3 

3 



Z 

« 


• 

(X 


u» 


m 



#- 



o 




Ui 




o 

O 



Ui 

• 


X 

44 


1-4 

« 

«• 




• 

3 





X 




o 

o 



z 

* 


1-4 

a 



in 

CU 



CX 

o 

O 

X 








<r 

<x 



o 

• 


a 

• • 


a 

a. 




Ui 

44 

mJ 

o 




in 








z 

• 



• o o 


►4 

Ui 




1- 

X 

x 




Ui 

w 




X 

X 

44 

9* 


« 


-j 

z z z 


CX 

i- 

* 



*4 

o 


• 

• 


UI 





Ui 

Ui 

X 

N 

X 

« 


<x 

3 ~~ 


<3 

in 






Ui 

o 

z 


X 

• • 




a 

a 


w 

X o 

• 


E 

OhH 

• 



> 



<x 

X 

X 

z 

o 

• 

o 

X • X 

• 



X 

X 

• 

• 

K *- 

• 


a 

XXX 

z 

z 

CL 

A 



1- 

o 

44 


44 

o 

Ui 

44 CX 44 1 

z 



o 

o 

z 

Z 

O 44 

• 


o 

p-» « <t 

o 

o 

Z 




<x 

u 

K 

X 

1- 

z 

a 

a 44 q 

o 





o 

O 

Z X 

« 


z 

OHH 

44 

44 


—J 



a 



Ui 

<r 


w 

Q 

44 



a 

o 

44 

44 


• 



in in 

H- 

K 

> 

O 




z 

UJ 

1- 

X 

X • 


<r <r 

1- 

• 


z 

z 

H- 

H* 

J X X 

♦ 


z 

Z Ui UJ 

o 

o 

(X 

IX 



-J 

o 


X 

Ui 

o * 

• 

h n h 

<x z 


ru 

ru 

<t 

<X 

w O 

« 


44 

P -4 <x (X 

Ui 

Ui 

Ui 

H 

* 


<x 

P -4 

X 


H- 

<x o 

> 

UiXN 

o o 



*• 

O 

o 

O * X 





in 

in 

3 

Z 

H- 


o 

H 

o 


44 

c <x 

X 


O 44 


ru 

^4 

o 

o 

44 X 



t- 

1 - cx cx 

z 


Ui 

o 

0 . 

* 

P -4 

o 

u. 

a 


•- 

<x 

z z z 

-J 



• 

• 

-J 

mJ 

K O 1- 

• 


z 

zoo 

o 

K 


o 


h- 

z 

X 


Ui 

A 

Z !- 

a 

44 44 44 




X 

X 



in 44 z 

• 


44 

« u. u n i-4 

<n 

(X 


i— a. 

o 

Ui 

in 

in 

O 

<x <X 

z 


X o 


Ui 

Ui 

1- 

H- 

Hh« 

• 

• 

o 

o 

a 

(X 

o 

<n 

cx 


o 

H 

•- 

3 

-J 

Ui 

3 

UJ Ui Ui 


o 

• 

A 

Q 

X 

X 

X <X O 

« 

<r 

GL 

x x x x a 

44 

u. 

44 

<x a 


44 

X 


o a x u. 

o 

rs> rsi m 

o j o tr 

X 

X 

<x 

Ui O X 

• 

Ui 


<X <X <X UJ 

U. 



H* 

X 

— J 





Ui K o 

A 

44 44 44 

x 


Ui o 

o 

K 

K 

1 - o 


X 

U 

UCCCL 


1- 

UJ 

U) Ui 

A 

u. 

X 

z 

u. in in 


in in in 

44 x in 



in 

in 

O -J X > 

<x 

O 

O -3 X 

o 

x 

3 

X 



•-4 

o 

o 

1- 

0 3 Ui 

X 


i- 


3 

H- 

(- 



<r 

• 



i- a 

o 

a. 

<r 

Ui Ui 

H- 





UJ _J 

Ui 

XXX 

X 1- 


in 

in 

Ui 

Ui 

X UI • 

• 

• 

• 

• a q <x 

• 

u. 

h- 

3 

* * 

44 

• 

• 

3 

• 

K UJ <3 

H- 

Ui Ui Ui 

<x in h 

94 



X 

<r cn u 

• 

u 

o 

O -J -J 

o 


3 

*• 

<x <x 

z 

o 

o 

Ui 

O O X z 

3 

H H H* 

h- <r o 

— 

— 

<x 

<x 

X O UJ 

• 

Ui 

z 

z o o * z 


O 

^4 

33 

44 

z 

z 

z 

z z u. <x 

o 

in in in 

in j z 

94 

a 

3 

3 

O Z X 

• 

X 


in 

A 



oox 



X 





x 

UI 







Ui 





H 1 

I- A 


X 

44 

<3 

<3 



X 


Ui 

UJ 

in 


X 

X X X z 



o 

in ui 

z 

Ui 

H- 

X 

X <t 

o 

<r <x 

4QUJ4 

cu 


X 

o 



x « <r ui z 

o 


X 

^ w 

o 

H 

X 

Ui 

UI Ui z u 

H 

h4 h 

<X Z UJ 3 

3 

<r 

<x 

JZ® <s 

<z 

z 

Z Z Z 1- X 

3 

X 

3 

3 3 

o 

44 

in 

1- 

i- -j in -j 

X 

Ui X N 

«- Ui 3 Z 

z 

3 

3 

X 1 - z z 

<r 


X “5 ^ * Z 

z 

z 

X 

X ^ 

z 

z 

z 

44 

44 ^ U. <t 

1 - 

AAA 

ru x in x 

X 

X 

csj 

O X X > 

<x 














A 

A 


A 

a 

« r- 

ID 














a a in 


• 

« 

• * cu 

00 


CD © ID MU 










• 

A 

♦ * c^ 


A 

94 

4 (J) (£ A 


o 

A 

in 

44 

aio 

A 

A 

cu 

cu 

4»TO» 

A 

44« 

A A • • 

• 

A 

• 

A CU * • 

• 

A 

rucutUH 

44 

A 


^4 m 

m 

m 



Ol • • 

r- 

• * • 

A A ~4 ~4 

A 

m 

A 

A t A A 

A 











A -4 <44 


A A cn 

<T CU 


94 





9 O O ® 9 Q 9 O 9 9.9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 
o-*tum'Tir>ior*aoo>9-«njn'Tir><DC'-oocn® — ajr)'«-inujr^aoOT®'-.rum'T_n 
HM^^M««^^^fl](u(unjfl)(unjnjn)njnnnnnnnnnnTTTTTT 


71 


H 

o 

03 

z 

M 


u 

0- H 

H 

x 

z > 

IM 

o 

O -J 

y A 

x 

*-» <x 

O- X > 


h r 

UJ CL 

X 

« <r 

CL X H 

o 

X •* 

cozou 

H * 

U Lu 

H O • X 

►H OC 

X ~ 

<Z ** ® O 

CL O 

UJ X 

EhhOLJ 

H 

o o 

O <X <X CD O 

a O 

UJ 

q: q j w 

o <X 

o a 

Ll. x 

u. u. 


H O 


0 c > 

CQIJOX 

X X 

o a 

h <r <i ♦ o 


HO O O <t H O ® H 
1 O • J ffi I Q </) O <X 
I z H O <x *-» UJUmQ. 
UJ<XUiOHXH(LO 
jno) H 3 oa a u. 
xaacia <iujo o 

• C UiZChhlL® 

It J^EhEU) w ♦ • 

luujoo** i lucrooo 
ddzzhao: 3 Q® z 

▼ 00 

o 


®HT 00 


ao m ▼ a> ao oo cu cu cu r- 

O CO ▼ LD r4 (D O) CO ® 


9000 t ▼ «T LD <-« LD 

m ld r- ob O) ® ® 


LO O') LO LO m 

ooooo 


® <0 

CD ® 

o a uj n 

CUJNQZ^OUJ 
OZCt:<X<XUUO ID 

JUJCOUi JQ.dK OO 
jCJDKKKroocn 
Cl OdZHtf Hhlt • ^ 

® LD 

cu ® 


H H H U. H H H 


▼ Cft 0J ID 

cn co ld n 


ID03C7> COOJ 

m ^ aj cu cu 


cu cn ao oo r* 

03 f* 00 «■« 00 


O^LDr^^OJ^nCU^LDCnLO^OLDLD 

® oi ao o- r- ld m ▼ m rn ^ i n lo 

«H I t 


LDr-ooa)®^cunxLDLor s -ooo>®^fun , TLDLDr^ooa)®^fucn <, rLDLDr-co 

▼ ▼^▼LOLOLDLDLDLDlDLDLDLOLDLOLDLDLOLOLDCDLOLDr-r-r-f^r-r^r^r^r^ 


72 


790- -90.72 112.20 J 

800- 571 

810- -90.72 112.20 I 

820- -99.15 110.39 V PATCH U 



co ® o ih u) h n ® 9in9T ▼ ^ ® m n cc t ^ 

^ co ao o <® o cu 'T 'rr-cuoof^ r*®Too ao <s ▼ ru r- r- p* ao cu 


m ® ® 90m nwoaniD tnioioao oo 10 a> a> ru cu to oj n 

<s ao in m cu in in r*^ o> ® ® q <© ® <s & ~* <** ^* -* ^oo>cd 

in 


o ® nj 


cn in ^ m id 

® at in ^ t p» «4 r* «4 <▼ o» cu cu o> ^ co ^ 

o* p* r*- n to ^ ▼ oo 'T 'r r* o> «* oi oi oo ao to ao p- nj m ao cu in 

*«**o>ino>ftjo>o>^nj , rm 4 r«<om'rcniDr-<UQDU>a)0>ai^ 
Ori^fuajr* o> o oo ao r- id u* ▼ t i m to to r- ao a> <» ©«^ai 

H ^ ^ 00 I I till »1 ^ ^ ^ ^ 

it it i f*i i i i i i 


99999900000000000000000000000000 

cn , riniDf s -cx30)®^run , TiniOf s -a30)®^cum , Tina)r^aoa)®^fucn , T 

030D0000«03»0>0i0i0>0)0>0>0>0>O®®®®®®®®®®^^^^^ 


73 


L 150- -125,54 30 . 98 

1160- -126.72 0.00 

L 170- 454.23 5 

L 180- 1 3 0 

L 190- 99.47 0.00 



SAMPLE OUTPUT 
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IV. OUTPUT DATA 


The program provides an output of. the flowfield after every NP march- 
ing steps (NP is user specified). A small portion of the output is presented 
here. It provides the flow variable information (density, Cp, <&, contravar- 
iant velocities, U, V and W, and the grid point location x, y and z) along the 
body surface (J * 2) and along the planes of symmetry (K * 2 and KMAX-1), 
Besides this printout for every NP marching steps, the code also prints the 
Courant number and maximun eigenvalues in the n and z crossflow directions, 
and the rms change in density. A detailed output of the entire flowfield in 
the crossplane (t),£) for all J and K can also be obtained at selected marching 
plane locations prescribed by the user in the input at card number 46 (see 
Section III). 

The code also stores the results for plotting purposes in Tape 1, and 
a separate plotting package is used to plot contours and other forms of visual 
output. 


V. RESULTS 

Four cases are presented to illustrate the capability of the code. 

1. Flow over an arrow-wing body at M,, * 2.96. a = 10.01°. 

2. Flow over a forebody configuration at = 2.5, 3 = 5° and at 
* 1 . 7 , a » 10 °, p * 5°. 

3. Flow over the entire Shuttle Orbiter geometry at M « 1.4, a * 

0 °. 

4. Flow over a realistic fighter configuration at different angles 
of attack and free-stream Mach numbers. 

For an attached shock at the nose (pointed nose configurations), a 
conical solution is generatd and used as an initial data plane for the non- 
conical marching calculation. If the nose is blunt and has a detached shock, 
then an unsteady full potential code 8 is used to generate the initial blunt 
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body solution (free-stream Mach number less than 1.5 is recommended). This 
initial flow calculation usually takes 40 to 60 iterations. At the beginning 
of each marching plane calculation, the* grid is generated using an elliptic 
grid solver. Usually, the grid solver requires 30 relaxation cycles to pro- 
vide an acceptable grid distribution. 

Figure 9 shows the pressure distribution on the surface of an arrow- 
wing configuration at location x/x » 0.8 for ■ 2.9, and <x * 10°. The im- 
provement In the prediction capability achieved using the wake treatment is 
illustrated. The dashed line represents the result from "no wake" treatment 
(assumes a flat plate behind the trailing edge) and the solid line represents 
the modifications to the pressure distribution once a zero jump in pressure 
across the wake cut is imposed. The solid line pressures on the body agree 
very well with experiments. 1 ^ Without a proper wake cut treatment, the over- 
all lift and drag forces and the pitching moment can be off by a considerable 
margin. 

Figure 10 presents the pressure distribution on a fully developed 
forebody 11 for ■ 2.5 and yaw angle S ■ 5°. Figure 11 shows the circumfer- 
ential pressure distribution for. the same body at x/i » 0.68 for = 1.7, a * 
10° and 0 * 5 8 . The experimental data 11 are also given in these two figures. 
The results show that the present prediction is in excellent agreement with 
the experimental data. 

Figures 12 to 22 give the geometry and the corresponding flow field 
solutions for the Shuttle Orbiter. The side view, cross-section and grid in 
the nose region are given in Fig. 12. The top view, cross-section, grid and 
chordwise cross-sections are presented in Fig. 13. It should be mentioned 
here that for most of the previous Space Shuttle calculations 1 ^’ 1 ^ the geom- 
etry has been modified by smoothing out the canopy and increasing the wing 
sweep angle from 45 to 55 degrees in order to avoid the subsonic bubble. 

Since the present method is valid for supersonic flow with embedded subsonic 
region, the realistic Shuttle Orbiter geometry was used without any 
modi fi cation . 

Figure 14 shows the surface pressure distribution along the leeward 
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plane of synmetry. At x s 170 in. which is the beginning of the canopy, the 
pressure Increases rapidly from C p 2 0.3 to 1.0, approximately. In the canopy 
region an MSR/TSR is formed and required three relaxation cycles to develop 
the solution. The circumferential surface pressure distribution at x * 240 
in. Is shown in Fig. 15. The experimental data is also given. in these figures 
and the agreement is very good. The surface pressure distribution along the 
wing leading edge is given in Fig. 16. It Is seen that the present calcula- 
tion agrees with the experimental data quite well along the entire wing lead- 
ing edge except in the vicinity of the wing-fuselage junction, where a vortex 
flow may exist. 

Figure 17 to Figure 21 present the Orbiter chordwise pressure distri- 
bution on the upper and lower wing surface at 7 * z/ bw » 0.471, 0.530, 0.641, 
0.780 and 0.887, respectively where bw Is a semispan defined in Fig, 12. The 
results show that the present predictions are in very good agreement with the 
experimental data, except in the region near the trailing edge of the upper 
surface at 7 * 0.471 and 0.53 span stations. Here again, a vortex flow or 
separation may be causing the discrepancies. Figure 22 shows the circumferen- 
tial pressure distribution for the Orbiter at x * 1120 in. It is noted that 
the pressure at the vertical tail and 0MS pods are well predicted. 

Figure 23 shows a supersonic fighter configuration with vertical 
tails. The free-stream Mach number, angle of attack and wing sweep angle for 
the calculation are summarized in Table 1. 

Figure 24 presents the surface pressure at various axial stations and 
the corresponding grid distribution for the wing body geometry. For this 
case, the MSR/TSR starts around x * 0.4 at the leading edge, and remains sub- 
sonic all the way to the end of the wing. Figure 25 and 26 show the circum- 
ferential pressure distribution in the vertical tail and wake region of the 
fighter-like configuration. The pressure on the vertical tail surface is 
given separately along the y-di recti on. The results clearly show that the 
present wake tretment provides the correct zero pressure jump condition across 
the wake. The chordwise pressure distributions from the center of the body to 
the tip of the wing are given in Fig. 27. Figure 28 presents the pressure 
contour on the upper and lower surface of the fighter configuration. 
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The lift and drag coefficients from the present calculations for this 
fighter model are also given in Table 1. The comparison with experimental 
data show excellent agreement. 


VI. CONCLUSIONS 

A nonlinear full potential method has been developed to analyze 
supersonic flows with embedded subsonic regions over some very complex con- 
figurations. A conservative switching scheme is employed to transition from 
the supersonic marching algorithm to a subsonic relaxation procedure. The 
predictions are in very good agreement with experiment data. 

Work is now progressing to simulate the multi body interaction between 
the Shuttle Orbiter and the external tank, and canard-wing fighter geometries. 
The present methodology will also be extended to treat all Mach number flows 
(fully subsonic flow, as well as. subsonic flow with pockets of supersonic 
region (transonic case)). 
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Fig. 3 Wake boundary condition. 
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Fig. 5 Geometry setup. 
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Fig, 12 Nose region geometry for 
Space Shuttle, 
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Fig. 14 Surface pressure distri- 
bution at leeward plane of 
symmetry. 
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Fig. 15 Surface pressure distri- 
bution around the body. 



Fig. 17 Shuttle Orbiter chordwise 

pressure distribution; M,, * 
1.4, a » 0.0°, 2 » 0.471. 



Fig. 18 Shuttle Orbiter chordwise 

pressure distribution; = 
1.4, a = 0.0° , 2 = 0.53. 



Fig. 16 Surface pressure distri- 
bution along the wing 
leading edge. 



Fig. 19 Shuttle Orbiter chordwise 

pressure distributi on ; M a * 
1.4, a = 0.0°, z = 0.641. 
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Fig. 23 Fighter-like configuration, 
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Table 1 Test cases for fighter-like configurations. 
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ABSTRACT 


Description of input/output data and operating 
instructions are presented for a recently developed 
supersonic full potential analysis program. Solution 
pre and post processors are also discussed for 
completeness. The latter includes a three dimensional 
finite difference boundary layer program interface 
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INTRODUCTION 


The supersonic full potential analysis method, 4 

in conjuction with three dimensional boundary layer analysis 
is currently being used to derive nonlinear supersonic cruise 
and maneuver designs to support aerodynamic advance configura- 
tion studies. The fundamental advance of full potential 
methodology relative to linear solutions is the ability to 
shock capture. The analysis consequently provides the means 
to achieve necessary conditions of shock weakening and 
separation elimination /management through appropriate re- 
design. In this regard, it provides a capability similar to 
relaxation solvers that are routinely used to develop 
efficient transonic flows. 

Supersonic full potential analysis also provides the 
ability to accurately assess the impact of sweep, thickness, 
and lift for a range of design space for which linear theory 
is unsatisfactory. 
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DISCUSSION 


The supersonic full potential (SFP) analysis code of 
Shankarl‘3 is applicable to arbit-rary wing-body-nacelle -tail 
arrangements from moderate supersonic Mach number to values of 
the hypersonic similarity parameter MiS'vO(l). The lower code 
limit is governed by the extent of embedded subsonic flow 
while the upper limit results from a breakdown in the isentropic 
assumption for strong shock waves. Also, since potential 
theory is irrotional, the modeling of vortices is not treated. 

The current version of the pilot code is operational on 
the CDC 875 and is stored as program library FPF7PL. 

CASE DESCRIPTION 


Four types of data are required to define a problem: 
a) header data describing mesh information, Mach number, angle 
of attack, aerodynamic coefficient reference quanities, center 
of gravity location, etc.; b) detailed geometric coordinates 
defining configuration cross plane contours ;'c) program update 
file directives defining code modifications, wake data, etc.; 
and d) job control directives defining program and input/output 
file allocations. 

GEOMETRY DATA 

A configuration is defined by several regions* of crossplane 
sections as indicated in figure 1. The number of patches 
(segments) defining a section is constant for a given region and 
typically increases from one region to the next. 

For the wing-body-vert ical case under discussion, a three 
(3) patch initial region, a six (6) patch center region, and 
and eight (8) patch final region as indicated on figure 2. Zero 
length patches are not permissable. Since the analysis is 
marching in nature, a complete geometry data set is not required 
to begin and partially process a problem. Appropriate use of 
restart solutions allows continuation of the analysis as new 
or modified geometry becomes available. 


*The overlap must be sufficient to encompass at least the final 
three (3) marching data planes of the prior region. 
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The format for a typical station is shown below. The 
group of cards is repeated for each station of a region. 

The last point of each patch (except for the last patch of 
a station) should have the same coordinates as the first 
point of the next patch. 

Card# Format Field Name Description 

A1 F15.6,I5 1 XI The x value (longitudinal) of 

this station. 

2 ISC1 The number of patches for this 

section. 1< ISC1 < 15 

The group of cards A2 thru A3 are repeated ISC1 times. 

A2 315 1 ITH Patch number < 15 

2 IPT1 Number of points in this patch. 

2 < IPT1 < 30 „ 

3 ND Mesh spacing parameter. Typically 

the same for all stations 'of a 
region 

The A3 card is repeated IPT1 times. 

A3 2F15.6 1 Y1 Vertical location of point 

(positive upwards) . Points 
start at top centerline. 

2 Z1 Spanwise location of point. 

Cubic spline interpolation is performed on input patch 
data to derive the boundary at the mesh points. Linear inter- 
polation is performed to define the boundary at a marching plane 
between input stations. 

Sample geometry data for the problem of figure 1 is 
presented on Table I and was developed using CDS^. 


*F0R SEGMENT AB 

0 Equal space 

1 Bunch near A 

2 Bunch near B 
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TABLE I GEOMETRY DATA 
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HEADER DATA 
* 

The present problem requires a minimum of four (4) 
marching solutions to complete the analysis of the three 
region sample geometry of figure 1 - one each for region 
one and two and two solutions for region three. The latter 
is a consequence of increasing the number of mesh points 
on the vertical tail as its local span increased with marching 
distance. It also illustrates the wake restart procedure. 

Each solution has a different set of header instructions 
for describing grid parameters, wake information if pertinent, 
restart directions, and number of mesh points for each patch 
of the region. This information precedes the geometric data 
discussed in the previous section. The header data used for 
the sample problem and a description of the various parameters 
is presented on Table II. 

The last header data set is coordinated with the wake 
update file described in the next section. The pertinent 
nomenclature is depicted below. 

KWAK* ** RESTART 

1 Standard (geometry) 

2 . Wake 


NU0* 



KWKED* 


NOTE: Reindex K of VORT (K) to allow for increase of 

points for patches 4 and 5 


* TAPE 5 

** Update file 
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TABLE II HEADER DATA 


Original fags is 

OF POOR QUALITY 


130= 

175 

NHAX 15 

1 10* 

15 

UMAX 

120= 

33 

KMAX 

138= 

2 

NRM 

148= 

13 

nub 

15^- 

2 

NP 

188= 

18 

KWKST 

170= 

36 

KWKED 

180= 

38 

NCCN 

198= 

38 

NITER 

290= 

S 

NSPTI 

210= 

l 

ITERGE 

229= 

S. 8 

CFLIN FI 8. 5 

230= 

2.5 

DZTAIN 

240= 

3.8 

DZMAX 

250= 

1.0 

DZMIN 

268= 

1.6 

FSM 

270= 

4.46 

ALFA 

280= 

68. 

THTO 

290= 

0.1 

DETA 

300= 

0.1 

DXI 

310= 

1.5 

DZTA 

320= 

15.00 

ZTAI 

338= 

371.3 

XEM> 

340= 

t. 

AMU1 

353= 

0. 

AMU2 

380= 

614. 

XVAKE 

370= 

117. 

ZVAKE 

380= 

1.0 

CHL 

360= 

85.3 

PTNOSE 

400= 

0 . 

YSHIFT 

418= 

580. 

XO 

420= 

a. 

YO 

438=125288. 

tWt 

448=234.627 

ALL 

453= 

1.75 

OMEGA 

480= 

T • 

PLA*€OS 

478= 

T 

NUGRID 

488= 

T 

I READ 

490= 

F 

RPLAI€ 

506= 

F 

TAPER 

518= 

T 

TAPE* 

528= 

F 

TAPE8W 

530= 

T 

FORCE 

540= 

18 

550= 

0. 

00.0 

588=180. 

150. 280. 

570= 

3 

ISC 

580= 

8 

10 15 

590= 

0 



NO. OF AXIAL STEP 

NO, OF PTS IN NORMAL 0tR.<28 

NO. OF PTS IN CIRCUM DIR: ISC! +ISC2-1 + I5C3+K81 

NO. GRID REGI0NS<7: B&WB-2, WBV-3, WBNV-5 

SECOf© SHARP EDGE K. 

OUTPUT FOR EVERY NP STEPS 
WAKE ST#?T K. 

I.MV C C\JT\ If 

noTcoc starting sol. steps. 

NO. OF GRID ITERATIONS. 

NO. OF ZTA FOR FLOW FIELD OUTPUT. 

NO. OF GL08 ITERATION. 

CFL NUMBER 

IF>0 FIXED STEP SIZE. IF<3 CFL NO. 

MAX. AXIAL STEP SIZE 
MIN. AXIAL STEP SIZE 
FREE STREAM MACH NO. 

ANGLE OF ATTACX-OEG. 

OUTER 8CLNDAPY-DEG. 

STEP SIZE IN ETA DIR. 

STEP SIZE IN XI DIR. 

FIRST STEP AFTER CONE START, SOL. 

STARTING ZTA >3*OZTAIN. 

END ZTAOVAX INPUT ZTA-OZTAIN. 

1 .‘FIRST ORDER, 2: 2ND CRC€R. 

8: FIRST ORDER, 1 :2N) ORDER. 

W <*E MINIMA ZTA. 

WAKE MINIMUM . 

GEOMETRY SCALE FACTOR 

AXIAL GEOMETRY SHIFT FOR ZTAD0. 

VERTICAL ISOMETRY SHIFT 
AXIAL C.G. ZTA. 

VERTICAL C.G. 

REFERENCE AREA. 

REFERENCE LENGTH. 

RELAXATION FACTOR. 

PLANE OF SYM,? 

SNERATE GRID? 

INPUT GEOETRY? 

R-MARCHING? 

RESTART DATA FROM TAPE? 

WRITE RESTART DATA CN UNIT 2*4? 

WRITE SUBSONIC RESTART DATA ON UNIT 8? ' 

CAL CUL ATE FORCES? 

GRID REGION TERMINAL K,5I5; ISC1+ISC2-I+. .+ISCN 
00.0 00.0 • POLAR tf£LE-0£S;SF13.4 

258. 308- 358. FLOW FIELD OUTPUT ZTA 

NO. GEOM. SEGfENTS 
NO. MESH PTSX SEGMENT 

NAJS 


a) Region 1 
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TABLE II CONTINUED 


100= 

125 

110= 

15 

120= 

57 

138= 

2 

140= 

13 

150= 

2 

168= 

28 

170= 

47 

180= 

30 

180= 

38 

208= 

6 

218= 

I 

220= 

5.9 

230= 

2.5 

248= 

3.0 

253= 

1.3 

260= 

t .6 

270= 

4.46 

280= 

68. 

290= 

0.1 

300= 

0.1 

313= 

1.5 

320= 

371.0 

338= 

611.0 

340= 

1 . 

358= 

0. 

360= 

614. 

378= 

117. 

380= 

1.8 


85.0 

400= 

0.0 

410= 

500. 

420= 

0. 

438=125288. 

440=234.627 

450= 

1.75 

460? 

T 

470= 

T 

488= 

T 

480= 

F 

580= 

T 

510= 

T 

520= 

F 

530= 

T 

548= 

32 

550= 

0.0 

560=400. 

570= 

6 

580= 

84 

590= 

0 


NO. OF AXIAL STEPS 

NESH PTS IN NORMAL DIR. <25 

MESH PTS IN CIROJM DIR: ISC1+ISC2-1+. .+ISCS+K81 
NO; GRID R£GIQNS<7; B4WB-2, WBV-3, WBNV-5 
SECOND SHARP EDGE K. 

OUTPUT FOR EVERY NP STEPS. 

WAKE START K; ISC1+ISC2-1. . . +ISC5 
WAKE END K;ISC1+ISC2-1+...+ISC7 
NOC0»«£ STARTING SOL, STEPS. 

NO. OF GRID ITERATIONS. 

m. OF ZTA FOR FLCW FIELD OUTPUT. 

NO. OF GLOB ITERATION. 

I^FdSd'sTE? S1ZE.IF<0 CFL NO~ 

MAX. AXIAL STEP SI£ 

MIN, AXIAL STEP SIZE 
FREE STREAM MACH NO. 

ANGLE OF ATTACX-OEG. 

OUTER 8CUN5ARY-OEG. 

STEP SIZE IN ETA DIR. 

STEP SIZE IN XI DIR. 

FIRST STEP AFTER COfE START. SOL. 

STARTING ZTA. ^ xTV1 

EM3 ZTAOtAX INPUT ZTA-DZTAIN. 

Is FIRST ORDER, 2: 2ND ORDER. 

0:FIRST ORDER, 1 s 2MC> ORDER. 

WAKE MINIMUM ZTA. 

WAKE MINIMUM Z. 

GEOMETRY SCALE FACTOR 

AX$E GEC^TRY SHIFT FOR ZTAD0. 

VERTICAL GEOMETRY SHIFT 
AXIAL C.G. ZTA. 

VE RTICA L C.G. 

REFERENCE AREA. 

reference length. 

RELAX I AT I CN FACTOR. 

PLAfC OF SYM.? 

GENERATE GRID? 

INPUT GEOMETRY? 

R-MARCHING? , 

RESTART DATA FROM TAPE? 

^ITE RESTART DATA ON UNIT 21 4? 

WRITE SUBSONIC RESTART OAT A CN UNIT 8? 
r^injLAT^ FORCES? 

SlDffiGION TERMINAL K,5I5j ISC1 + ISC2H+. .+TSCN 
an 9 00.0 - PCLA* ANa£:Sr10.4 

553. 575 600. FLOW FIELD OJTPUT ZTA 

w NO. GEOM. SEGMENTS _ 

S 1 1 NO. «SH PTS/SEGMENT 


IS 11 
NAUS 


b) Region 2 


ORIGINAL PAGE IS 

OF POOR QUALITY 
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TABLE II CONTINUED 


ORIGINAL PAGE IS 
OF POOR QUALITY 


108= 

50 

IFAX 

110= 

15 

UMAX 

120= 

63 

KMAX 

1 30= 

3 

fFM 

140= 

21 

NU0 

150= 

2 

IF 

168= 

24 

KWKST 

179= 

52 

KWKED 

180= 

30 

NCCN 

160= 

30 

NITER 

208= 

3 

NSPTI 

210= 

1 

ITERS 

220= 

5.0 

CFL IN 

238= 

2.5 

DZTAIf 

240= 

3.0 

DZMAX 

2S3= 

1.0 

DZMIN 

268= 

1.6 

FSM 


270s 

280= 

280 s 

380= 

318* 


4.46 

60. 

9.1 

9.1 

1.5 


320= 611.9 
338= 685.9 
1. 


340 

350s 

360s 

378= 


300s 

490s 

410= 

429= 


0. 
614. 
117. 

1.9 
95.9 
0.8 
500. 
0. 


438> 125280. 
449=234.627 


450= 

460= 

470= 

480= 

460= 

500= 

5!0= 

520= 

538= 

549= 

950= 


1.75 

T 

T 

T 

F 

T 

T 

F 

T 

21 

70.9 


560=625. 
578= 3 

580= 04 

560= 0 


ALFA 

THTO 

DETA 

DXI 

DZTA 

ZTA1 

XEND 

AMUt 

AMU2 

XWAKE 

ZVAKE 

CH. 

PTNGSE 

YSHIFT 

XO 

ro 

W* 

ALL 

QtCGA 

PLANEQS 

NUGRID 

I READ 

RP U*€ 

TAPER 

TAPEW 

TAPE8V 

FORCE. 

38 80 00 

0.00 

658. 675. 

ISC 

85 18 84 


15 NO. OF AXIAL STEPS 

f£SH PTS IN NORMAL DIR.<2S 

«SH PTS IN CIRCUM DIR? ISC1+ISC2-H*. .+ISC8+K81 

NO. GRID REG I QNS<7 ; 8&WB-2 . wev-3, V8NV-5 

TAIL EDGE K; ISC1+ISC2-1 +. . . +ISC4 

OUTPUT FOR EVERY NP STEPS 

WAKE START K;ISC1+ISC2-1+...+ISCS 

WAKE END K: ISC1 +TSC2-I *. . .+ISC7 

NO. CONE STARTING SOL. STEPS. 

NO. OF GRID ITERATIONS. 

NO. OF ZTA FOR FLOW FIELD OUTPUT. 

„ NO. OF GLOB ITERATION, 

F18.S CFL NUMBER. 

IF>0 FIXED STEP SIZE. IF<2 CFL NO. 

MAX. AXIAL STEP SIZE 
MIN. AXIAL STEP SIZE 
FREE STREAM MACH NO. 

ANGLE OF ATTACX-CES. 

OUTER BOlf'CARY-DES . 

STEP SIZE IN ETA DIR. 

STEP SIZE IN XI DIR. 

FIRST STEP AFTER CX» € START. SCL. 

STARTING ZTA. 

END ZTAOIAX INPUT ZTA-OZTAIN. 

1: FIRST ORDER, 2: 2ND ORDER. 

0:FIRST ORDER, 1 s2ND ORDER. 

WAKE MINIMUM £TA. 

WAKE MINIMUM Z. 

GEOMETRY SCALE FACTOR 

AXIAL GECMRTRY SHIFT FOR ZTAD0. 

VERTICAL GEOMETRY SHIFT. 

AXIAL C.G. ZTA. 

VERTICAL C.G. 

REFERENCE AREA 
REFERENCE LENGTH. 

RELAXATION FACTOR. 

PLANE OF SYH,? 

GENERATE GRID? 

IN=UT GEOMETRY? 

R-MARCHING? 

RESTART DATA FRQ1 TAPE? 

WRITE RESTART DATA ON UNIT 2 i 4? 

WRITE SU8SCNIC RESTART OAT A ON UNIT 8? 

CALCULATE FORCES? 

GRID RESIGN TERMINAL 
80.8 00.0 “ 


04 IS 
NAJS 


15 11 


K,5IS; ISC1+ISC2-1+.. + ISCN 
PQ-Afl AMa£}^18.4 
FLOW FIELD OUTPUT ZTA 
NO. GECM. SEGJ-ENTS 
NO. fESH PTS/Sc GHENT 


c) Region 3, Solution 1 
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TABLE II COMPLETED 


CniGiNAL PAGE IS 

OF POOR QUALITY 


100= 

1 10 * 

120 * 

130* 

140* 

150= 

160* 

170* 

188= 

168* 

290* 

210 * 

228* 

230= 

240= 

253= 

268= 

273= 

288= 

268 = 

380= 

310= 

328= 685.0 
338= 884.0 
348= 

358= 

380= 

370= 

388= 

360= 

408= 

418= 


75 NMAX 15 

15 JMAX 

69 KHAX 

3 NRM 

24 NU0 

2 NP 

38 KWKST 

53 OKED 

30 MCON 

38 NITER 

5 NSPTI 

1 I 

5.8 CFLIN F10.5 

2.5 DZTAIN 

3.0 Dzmx 

1.0 DZMIN 

1.5 FSM 

4.46 ALFA 

68. THTO 

0. 1 DETA 

0. 1 DXI 

1.5 DZTA 

5.0 ztai 

J4.0 XEM5 

t. AMJt 

0. AHU2 

614. XWAKE 

117. ZWAKE 

1.0 CH_ 

85.0 PTNOSE 

0.0 YSHIFT 

509. XO 

0. YO 


430=125288. 



440=234.627 


ALL 

450= 

1.75 


CJEGA 

460= 

T 


PLANEOS 

470= 

T 


NUGRID 

488= 

T 


IREAD 

460= 

F- 


RPLA4E 

560= 

T 


TAPER 

518= 

T 


TAPEW 

520* 

F 


TAPE8W 

538= 

T 


FORCE 

540* 

24 

44 

98 a 

550= 

78.0 


0.00 

560*708. 

725. 750 

578= 

8 


ISC 

588* 

04 

05 

18 

590= 

1 

04 

04 a 


NO OF STEP 

MESH PTS IN .NORMAL DIR. <26 

MESH PTS IN CIRCUM DIR: ISC1+ISC2-1 +. . +ISC8+K31 

NO. GRID REGICNS<7 ; B&W8-2, W8V-3, W6NV-5 

TAIL EEKS K; ISC 1 +ISC2-1 +. . . ♦ISC4 

OUTPUT FOR EVERY NP STEPS 

WAKE START K: ISC1+ISC2-1+. . • +ISC5 

WAKE EM) K; ISC1+ISC2-1+. iA +ISC7 

NO. 0>£ STARTING SO.. STEPS. 

NO. GRID OF ITERATIONS. 

NO. OF ZTA FOR FLOW FIELD OUTPUT. 

NO. CF GLOB ITERATION. 

IF>0 N FI^D*STEP SIZE. IF<0 CFL NO. 

MAX. AXIAL STEP SIZE 
MIN. AXIAL STEP SIZE 
FREE STREAM MACH NO. 
fiNGLE OF ATTACK-OEG. 

OUTER BOUNDARY -DEG. 

STEP SIZE IN ETA DIR. 

STEP SIZE IN XI DIR. ^ 

FIRST STEP AFTER C1>£ START . SOL. 

STARTING ZTA. . _ . . 

EM) ZTA<MAX IMVT ZT A-QZT AIN. 

1 -.FIRST CRD£R,2:2M> GRDER. 

0: FIRST ORDER, 1 :2ND ORDER. 

WAKE MINIMUM ZTA. 

WAKE MINIMUM Z. 

GECftTRY SCALE FACTOR 

AXIAL GEOMETRY SHIFT FOR ZTA>8. 

VERTICAL GEOMETRY SHIFT. 

AXIAL C.G. ZTA. 

VERTICAL C.G. 

REFERENCE AREA 
REFERENCE LENGTH. 

RELAXATION FACTOR. 

PLAICE CF SYM.? 

GENERATE GRID? 

INPUT GEOMETRY? 

R-MARCHING? _ A ___ 

RESTART DATA FROM TAPE? t 

m Sit st 

. *iscn 

®- 9 ‘ 

7751 HO. SEEM. SB?£HTS_ 

7 15 IS 11 NO. MESH PTS/SEGMENT 

04NAJS, NPTCHA, NOCPTA, M>TCH8, NCCPTB 


d) Region 3, Solution 2 
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UPDATE FILE DATA 


The purpose of this data is to define specialized 
information logic which has not been incorporated into the 
pilot code. It is anticipated that this type of input will 
be progressively eliminated as a production status evolves. 

Two basic update files (see table III) are typically used. 
One for a standard (i.e. geometry) restart and a second for a 
wake restart. The former defines certain potential averaging, 
leading edge mesh, nacelle directives, etc. In addition, the 
wake restart file describes the vorticity distribution which 
is output by the prior (for the present case the third) 
solution. 

The update files also contain modifications to the pilot 
code which have not been incorporated permanently. These 
directives do not need to be changed (but must be included) 
for a new problem. 


105 



TABLE III UPDATE FILE DIRECTIVES 


1 00** I DENT HTOWIR6 
I 10**1 MAIN. 14 

120*0 taped.- output pressure distributions at specified z 

130**0 MAIN.SS.MAIN.72 

!sl=C RESTART DIR£CT[vE;KWAK=1-SCLID BOUNDARY. -2- WAKE 
IS0-C 
1 70* 

180**0 
1 90*C 
200 - 
210 


<WAK* 1 

INVSETA. 1 7S. INVSETA. 1 7S 

AVERAGE PHI AT MAX Z „ 

IF ( ZTAI .LT. 371.0 ) PHK I ,2.KTE«»>* 
*(PHI« 1,2.KTEMP-H l+PHII I.2.KTEMP.-' >>/2-0 


220**0 NFORCE. 1 SI 


18 FORMAT ( 2X. 'PITCH MOHENT* ' . £12. 4. 3X. 'CM-\EI2.4, 
INPUT. 27 

DO 10 NAA-I.S0 


230 ; 

240**0 
263= 

2S0*«O METRIC. SI .METRIC. S2 


270-C 
280-C 
290*«ECF 


CONST. FOR. ZTAI SHOULD BE CHANCE AT VERY THIN. SHARP 
LEADING EDGE USING fCTRIC So- 79 


a) Region 1 


I 00-* I CENT HROWIRE 
I 10**1 MAIN. 14 

1 20=C TAPEOs OUTPUT PRESSURE DISTRIBUTIONS AT SPECIFIED Z 

130**0 MAIN. SS. MAIN. 72 
1 40 *C 

1S0-C RESTART DIRECTIVE: KWAK*1 -SOLID BOUNDARY , -2-WAKE 
I60=C 

1 70* KWAK-1 

180**0 NFORCE. 1 SI 

190= 18 FORMAT (2X. PtTCH MOMENT*' , £1 2. 4. 3X. 'CM*'.£12.4, 

200**0 GEOM. I7S 

210*C INCREASE NO. CF ,“ESH POINTS IN WING L.E. REGION 
220* DS-SSS/(NPT( IS)*3.SJ 

230**0 GEOM. 1 86 

240=C INCREASE NO. OF MESH POINTS IN WING L.E. REGION 
2S0* DS=SSS/(A0*3.S) 

2S0=*O INPUT. 27 

270* DO 10 NAA-I.S0 

280**0 METRIC. SI . METRIC. S2 

290 =C CONST. FOR ZTAI SHOULD BE CHANGE AT VERY THIN. SHARP 

330-C LEADING EDGE USING fETRIC So- 79 
31 0**£OF 


b) Region 2 


I 00** l DENT HROWIRE 

{ ^l=C 1 MAIN * TAPED, OUTPUT PRESSURE DISTRIBUTIONS AT SPECIFIED Z 
130**0 MAIN. 5S. MAIN. 72 
t 40=C-«* 

1S0*C*» RESTART DIRECT I VE;<WAK-1 -solid BOUNDARY , =2-WAK£ 
i63=c** 

179* KWAK=1 

180**0 MAIN. 300 

190=0 DEFINE WING T. £. EQUATION 
200* ZWAKE ! * \ . S622*ZTAl - 1 026. S3 

21 0=*O HA I N . 3 1 5 , MA I N . 3 1 6 


229= C 
230* 
240= 
2S0=*O 
2S0=C 
2 / 0 = 
280=*D 
2S0=C 
300= 
.310* 
320= 
339= 
340=*O 
3S0= 
360 -*0 
370=C 
380= 
390 =*0 
400=C 
410 = 
422 = *0 
430 = 
440 - *0 
4§0=C 
463 = L 


SPECIFY HAXIMJM W**E Z<I.£* Wl M3 HALF SPAN) 
IF (ZWAK.El .GT. 306.53) XVKST I =X7EItP 
“ 306.53) GO TO 54 


ORIGINAL PAGE IS 

OF POOR QUALITY 


FOR WAKE RESTART 
GO TO 138 


4/3**£GF 


IF ( Z WAKE I .GT. 

MAIN. 3S0. MAIN. 361 
OUTPUT VCRTICITY 
IF <NI .NE. 291 
INVSETA. 1 75. INVSETA. I 76 
AVERAGE PHI AT OR NEAR WING TIP 
IF t ZTAI .GT. S4S.0 .ANO. ZTAI .Lt. 717.0) PHI ( 1 , 2. :<TEW» > = 

*( PHI t 1 ,2.KTErt»*l )*PHI ( 1 .2.KTEMP-I ) >/2.0 
IF( ZTAI .GT. S4S.0 .AND. ZTAI .LT. 717.0) PHI < 1 . 2. KTEMP*1 )> 
*( PHK 1 , 2. KTEMP ) *PH I ( 1 ,2.KTEMP*2) )/2.3 
NFORCE. I SI 

18 FORMAT <2 X. PITCH mOMCNT* ' 

GEOM. i ;s 

INCREASE NO. OF MESH POINTS IN WING L.E. REGION 
DS-SSS/I NPT ( IS ) *3.5 > 

G£CM. i 86 

INCREASE NO. OF MESH POINTS IN WING L.E. REGION 
DS-SSS/ ( A0*3. S ) 

INPUT. 27 

DC 10 NAA-I S0 
METRIC. SI METRIC. S2 

CONS' . FOR ZTAI SHOULD BE CHANGE AT «ERY THIN. SHARP 
LEADING EDGE USING METRIC SS - ’9 


.E12.4.3X. "CH= ',£12.4, 


c) Region 3, Solution 1 
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ORIGINAL PAGE IS 

OF POOR QUALITY 

TABLE III COMPLETED 


!00=*ID£NT WS1WIRE 
118=*I MAIN. 14 

120=0 TAPES: OUTPUT PRESSURE DISTRIBUTIONS AT SPECIFIED Z 

138=*0 MAIN. 55, MAIN. 72 

140=0 

150=0 RESTART DIR£CTIVE;KVAK=1-SCLID BOUNDARY, =2-WAKE 
168=0 KVED1-MAX U.S. WAKE K; KWKST1-MIN L.S. WAKE K 
170=0 D783-15A 55 DEG W.T. MODEL DESIGN; M=1 . S,ALPHA=4. 48 

183=0 USE NACELLE OFF FILE 0S5C, ITER=2S WITH REINCEXED K 

130=0 

200= KWAK=2 

210= KVED1=3S 

220= KV*ST1=S2 

230= VORT(30)=28.5t 

240= YORK 31 ) =28. 27 

250= VCRT(32)=27.S2 

250= Y0RT(33)=27. 49 

279= VORT( 34) =28.03 

280= Y0RT(3S)=23.31 

290= VCRT(36)=19, 83 

380= V0RT<52)=-19.84 

310= VCRT<53)=-23. 31 

320= V0RT(S4)=-2S, 04 

338= V0RT(55)=-27.49 

343= V0RT(56)=-27.52 

350= V0RT(57)=-28. 57 

360= VGRT ( 58 )=-28. 82 

370=40 MAIN. 388 

380=0 DEFINE WING T.E. EQUATION 

398= ZVAK£1*I . 8622*ZTA1 -102S.53 

400=*O MAIN. 31 5, MAIN. 316 

410=0 SPECIFY MAXIMUM WAKE Z< I.E. WING HALF SPAN) 

428= IF (ZVAKE1 .GT. 306.53) KWKST 1 =KTEMP 

438= IF CZWAKE1 .GT. 306.53) GO TO 54 

443= *0 MAIN. 360, MAIN. 381 

450=0 OUTPUT VORTICITY FOR WAKE RESTART 

468= IF (NI .NE. 48) GO TO 138 

470=*O INVSETA, 175, INVSETA. 178 

488=0 AVERAGE PHI AT OR NEAR WING TIP 

490= IFIZTAl .GT. 648.0 .A*©. ZTA1 .LT. 717.0) PHK 1 ,2,KTEMF)= 
500= *<PHI( 1 ,2,KTEMP+1 )+PHI< 1,2,KTEMP-1 ) )/2.0 

513= IF1ZTA1 .GT. 648.0 .AND. ZTA1 .LT. 717.0) PHK 1 , 2,KTE?F-M )* 
520= 4<PHI(1,2,KTEW > )+PHI(1,2,KTE>f>+2))/2.0 

530=0 AVERAGE PHI AT OR !«AR VERTICAL TIP 

540= IF1ZTA1 .GT. 760. > PHK 1 ,2,24)=<PHK 1 ,2,23)+PHK 1,2,25) )/2. 

550= IF1ZTA1 .GT. 768. ) PHK 1 ,2.25 )=PHI( 1 , 2,24) 

560=40 M^CRCE.ISI 

570= 18 FORMAT (2X, 'PITCH MOMENT- '.El 2. 4, 3X, 'CM=',E12.4, 

588=40 GRID. 133 

598=C INCREASE GRID CELL AREA IN VERTICAL/WAKE CORNER 

ff&rGRiSn^ 75 ' 

620=0 CONVERGE- MORE SLOWLY TO AVOID GRID DIVERGANCE 

638= IF<W .EQ. I ) OP =.005 

648= IF1M4 .GE. 2) OP*. 01 

650=40 INPUT. 27 

688= 00 10 NAA= 1 , 50 

670=40 METRIC. 51, METRIC. 52 _ _ 

688=0 CONST. FOR ZTA1 SHOULD BE CHANGE AT VERY THIN, SHARP 

690=0 LEAOING EDGE USING METRIC 56-79 

788=*E0F 


d) Region 3, Solution 2 
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JOB CONTROL 


The final type of information required to process a 
problem is program/update and input/output file declarations 
These directives change from solution to solution in order 
to properly process and save data. 


The pertinent submit files for the CDC 87 5 under N'OS 
Operating System 2.1 are presented on table IV for the sampl 
problem. The file functions are: 


TAPE 1 
TAPE 2 
TAPE 3 
TAPE 4 
TAPE 5 
TAPE 7 
TAPE 8 
TAPE 9 


Last X Plot Data 

Output Marching Plane Restart Data 
Input Marching Plane Restart Data 
Output Marching Plane Restart Data Backup 
Input Data 

Input Subsonic Region Data: 30 Step Limit 

Output Subsonic Region Data: 30 Step Limit 

Output Pressure Distribution at Specified Z 
for Post Processor 



ORIG^AL PAGE !S 

TABLE IV SAMPLE JOB CONTROL DIRECTIVES OF POOR QUALITY 


i 08=SS5A, £C700 , T300, P3. 

1 1 0=USER, D8236, XXXXXX. 

1 20=CHARGE, *011, XXXX. 

1 30=ACCT ( 8CNNER ST 1 8422401 *01 1 XXXX ) 

140=RFL,EC=708. 

I S3=ATTACH, 0LDPL=SFP7PL/UN=0883S. 
160=GET, INP=US5A. 

173=UPDATE, I=INP. 

1 80=FTNS, I , OPT =2, LCM*I , L=Q, PL= 20000. 

1 90=ATTACH, TAPES=AT55A. 
200=PURGE,PLS5A/NA. 
210=OEFrNE.TAPE1=PLS5A, 
220=PL)RGE,RS5Al/NA. 

230=OEF [NE, TAPE2=R55A1 . 

240=PURGE, RS5A2/NA. 

250=DEFINE, TAPE4=R55A2. 
260=PURGE,PS5/NA. 

270=OEFINE, TAPE9=PS5. 

280 =A TT ACH , L IB 1 =SFP7LGO AN=O083S . 
230=CQPYL,LIB1 ,LG0,LG01 . 

300=LGO1 . 

3 1 0=*EOF 


1 00=S558 , EC700 , T500 , P 3. 

1 1 0=USER, D0236, XXXXXX. 

120=CHARGE, *01 1 , XXXX. 

1^*ACCT(8CNNER ST1 8422401 *01 1 XXXX) 

\ 5^TTACH, O^L=SFP7PL/UN=O083S . 
1S8=G£T,IW 3 =U558. 

cptS; lcm=i , 1=0, ^-=20000. 

1 Q0=ATTACH , TAP£5=ATS5B . 

200=ATTACH, TAPE3=RS5A1 . 

2 1 3=PURG£ < PLS5B/NA * 

228=0EFINE, TAPE! =PL55B. 

230=PUROE, R55B1 /NA. 

248=0EFINE, TAPE2=RSSB1 . 

2S0=PURG£, RS582/NA. 

268= DEFINE, TAPE4=RSSB2. 

278= ATTACH, TAPE9=P55/H=A. 

288=SKIPEI, TAPES. 

298=ATTACH, LIBt =SrP7L0O/UN=D083S. 
308=COPYL, LIB1 , LGO, L601 . 

310=0331. 




a) Region 1 


b) Region 2 


1 80=SSSC, EC700, T200, P3. 

1 1 3=USER , D0236 , XXXXXX. 

1 28=CHARGE, *81 1 , XXXX. 

130=ACCT(BOf*£R ST1 8422401 *01 1 XXXX) 
140=W r L,EC=700, 

1 50=ATTACH, OJ3PL=SFP7PLAJN=C883S. 
168=GET, I!*»=US5C. 

173=CPDATE, I=INP. 

188=FTN5, 1 , OPT =2, LCJ4= l , L=0, PL=28880. 
1O0=ATTACH, TAPS5=AT55C. 

200=ATTACH, Ttf»E3=RS581 . 
210=PURGE,PL5SC/NA. 

220=DEFINE, TAPEJ=PLS5C. 
230=PURGE,RS5C1/NA. 

240=OEFINE, TAPE2*RSSC1 . 
250=PURG£,RS5C2/NA. 

268=C€FIh€, TAPE4=RS5C2. 

278=ATTACH, TAPE3=PS5/H=A. 
280=SKIPEI,TAPE9« 
290=ATTACH,LIB1=SFP7L0OAJN=O083S. 
308=COPYL,LIB1 ,LGO,LGOt . 

310=0301. 

320=*EOF 


c) Region 3, Solution 1 


1 00=S55D , EC700 , T300, P3. 

1 10=USER,D0236, XXXXXX. 

120=CHARG£,*31 1 , XXXX. 

1 30=ACCT ( BONNER ST1 8422401 *01 1 XXXX > 

1 40=RFL , EC= 708 . 

1 53=ATTACH, OLDPL=SFP7PL/UN=O0835. 
168=G£T, INPHJS5D. 

170=tPDATE, 1= INP. 

I80=FTNS,I,QPT=2,LCM=l,L*O,PL=22083. 

1 98= ATTACH, TAPE5=ATS5D. 

208= ATTACH, TAPE3=R55Ct . 
210=PURSE,PL55D/NA. 

220=OEFINE. TAPE1 =PLS5D. 

230=PURGE, RESD1 /NA. 

240=C£FINE, TAPE2=RSSD1 . 

250=PURGE , RS5D2/NA . 

2S0=OEF Ih€ , TAPE4=RS5D2 . 

270=ATTACH, TAPE9=PS5/TI=A. 

*>oo — cvtoct TAPP9 

I^ATTACHi L IBT=SFP7LSO/LN=D083S. 
300=COPYL, LIB1 , LGO, LG01 . 

318=LGG1 . 

320=*EGF 


d) Region 3, Solution 2 
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ANALYSIS 


It is recommended that a general file nomenclature be 
developed prior to the analysis in order to achieve a con- 
sistent function/code/problem seven character descriptor. 
The following is provided as an example. 

FUNCTION FILE 


Submit 


Function 

\ 


Code Solutions 
/ 


SXYYYY A, B, C, D 




Case 


Input I 
Sectional Plot X 
Restart R 
Printed Output 0 
Pressure Plot P 


A review of the submit, header input, and update files 
should be completed before processing each solution. A 
check list which has been found to be useful in this regard 
is 


SOLUTION 

TAPE 5 

UPDATES 

Region 1 

Taper=F 

NAJS=0 

Modify as appropriate 

Region 2 

Taper=T 


Region 3 
Solution 1 

NRM-3, NU0 
KWST , KWKED , 
XWAKE 

Main. 360, VORT (k) Output 

Region 3 
Solution 2 

NAJS-1, etc 

KWAK* 2, Re index 
KWKED 1,KWKST1, and 


VORT (k) 
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Typical problem (one Mach number and angle of attack) run 
time on the- CDC 875 under OPT = 2 compilation is 


SOLUTION 

ZTA1 

J 

K • 

ITER 

CPU-SEC 

SEC/STEP 

Region 1 

0-370 

15 

33 

172 

121 

.7 

Region 2 

370-455 

15 

57 

34 

60 

1.76 

Region 3 

455-610 

15 

68 

62 

111 

1.78 

Region 4 610-685 

Solution 1 

15 

74 

30 

79 

2.65 

Region 4 685-8.02.5 

Solution 2 

..15 

80 

47 

133 

S04 

2 . 82 


111 


PRE- PROCESSOR 


i 






A pre-processor is available for evaluating input geometry 
and grid quality. It may be accessed. and executed using 

GET , SFPPLTX/UN=D0235 
SFPPLTX,XPFN , DOXXX 


The file XPFN is generated for ZTA1 by the supersonic full 
potential analysis as tape 1 when NMAX=0, TAPER=F. 

Typical grid program prompt response is 


. S085E*03 20 


<ret, 9fopitx 

,, '”!is08E?l^’ 4 JKee*9' 

[PLOT.^.iS ^ C 2 NTCURS 

- 3 "09 PRESSURE CONTOURS 

= 10 FOR COMPUTATIONAL GRID 

• 12 RQR maCH NO. PROFILES 

13 FOR PRESSURE PROFILES 


- 0 - 


- 8 - 


- 0 - 


0 


. 94091 86SE*03 
. 88S99023E-03 


•> 10 

XM!N= 0. XMAX= 

vmin= -.9061 9733E+03YMAX= 

XMIN, XMAX. YM!N, YMAX. . . 4F10.4 
< IS T HE PLOT IS X COORDINATE 
IS THE plot is y COORDINATE 
IS THUS M ACH NO. FOR iPLCT=12 
IS THUS PRESSURE FOR iPL0T=j3 
is OTHERWISE the PHYSICAL Y COORDINATE 
IS always the PHYSICAL X COORD I NATE 
vmtvjs q XMAX= .94091 865c 03 

v«IN= -'. 90619733E*03YHAX= .38S99023E*03 


XMIN= 0. 


- 0 - 


u 

YMIN= 

? 0 . 


XHAX= 


. 9061 9733E*03YMAX= 
400. -200. 


-0 

. 9409 1 86SE-03 
. 885990236*03 
200 . 


SI 


ORIGINAL PAGE SS 

OF POOR QUALITY 


Typical ourput is shown below. The grid is truncated 
at Z=400 and Y=±200 as specified in the input above 



ORIGINAL PAGE IS 
OF POOR QUALITY 


OUTPUT DATA 


Sample printed output data is presented on table V for 
region 3, solution 2. Standard tabulated data is produced 
every NP marching steps as defined in the header data. 

More detailed physical plane data can be output at 
specified stations using the parameter NSPTI and the 
tabulated stations following the true/false header data 
input. Cartesion coordinates, velocities, and mesh indices 
are indicated in the following sketch. The axial velocity 
component, U, is positive out of the plane of the paper. 


f Y,V 
k=2 



Crossplane data and tabulated surface pressure 
coefficient data at constant span stations (i.e. z) are 
output via tapes 1 and 9 respectively. These files are 
displayed using the pre and post processors described in this 
document. 
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POST PROCESSOR 

A graphics geometry and surface pressure coefficient 
output processor for the Techtronix 4114 has been recently 
developed to support advance design effort using the supersonic 
full potential analysis. It may accessed and executed using 

GET, SFPPLT/UN=D0235 

SFPPLT , FP IN , FPOUT , DOXXX 


The pertinent files are 

NAME TYPE DESCRIPTION 

SFPPLT ID Proc for Post Processor 

FPIN DA SFP Output Data Set < 5 Characters. 

UN=DOXXX 

FPOUT DA Created to Store Post Processor 

Output for Debug if necessary. 

The input file FPIN is generated by supersonic full potential 
analysis as tape 9 output. This file must be edited as follows: 

1. Delete *EOR's and repetitious results. 

2. Delete cone starting solution off front end. 

3. Add the following data at the beginning. 


(a) Title Card 

(b) Mach Number and Angle of Attack. One card, 
format 2E12 . 5. 

(c) Number of trailing edge points. One card, 
format 13. 

(d) Planform trailing edge data. Input trailing 
edge X from root to tip, for each section 
stored in FPIN. Format SE12.5 

4. Check end of data set to make sure the last x-station 
is fully represented. If ic is not, (examine previous 
x-station to see a fully represented example)' delete 
these incomplete data. 

5. If you do not want a particular z -value (3rd column) to 
be included in the data you will process, delete any 
reference to such data out of the 3rd column at the 
very end of the data set. 

Sample graphics post processor input data is presented on 
Table VI . 
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The following nine post processor options are available 
to the analyst. 


0 TO END 

1 TO DISPLAY PLANFORM 

2 TO DISPLAY 3-D VIEW OF WING Z/C 

3 TO DISPLAY 3-D VIEW OF WING CP S 

4 TO DISPLAY SECTION CP'S 

5 TO DISPLAY SECTION Z/C 

6 TO DISPLAY SPANWISE CL,CD,XCP 

7 TO DISPLAY PLANFORM ISOBARS 

8 TO CREATE BOUNDARY LAYER DATA FILE 


4 

The last creates a boundary layer analysis input file 
FPINBL for the upper and lower surface. This data must 
be edited to remove the first STOP card, update case 
identification, static pressure and temperature, and 
specified transition point location. 
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Patch 1 
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•igure 2. Region 1 Patching 
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Abstract 


An implicit, conservative treatment for the 
unsteady full potential equation in two -dime ns ions 
is presented. The method employs a local time 
linearization for density, and introduces flux 
biasing concepts based on sonic conditions for the 
generation of artificial viscosity to capture 
shocks without any overshoots. The boundary con- 
dition is treated implicitely using a splitting 
procedure consistent with the approximate factor- 
ization scheme. This allows for extremely large 
Courant numbers, even for nonorthogonal grid at 
the body. The method has application not only to 
unsteady problems, but also to generate the start- 
ing blunt body solution for a supersonic full po- 
tential marching code. Results are presented for 
flows over cylinders, spheres and airfoils. Com- 
parisons are made with available Euler and full 
potential results, and are in excellent agreement. 

I. Introduction 


Nonlinear aerodynamic prediction methods based 
on the steady form of the full potential equation 
are used ragularly for treating transonic 1 * 2 and 
supersonic 3 ” 6 flows over realistic wing-body con- 
figurations. Numerical algorithms to compute the 
unsteady form of the full potential equation is 
still in a developmental stage, and several re- 
searchers 7 ” 11 have recently made significant pro- 
gress in this area. There are several issues 
involved in the construction of a robust and effi- 
cient numerical algorithm for the unsteady full 
potential equation. They are: 1) proper treat- 

ment of boundary conditions in a nonorthogonal 
grid system, 2) correct formulation for the pro- 
duction of artificial viscosity to capture sharp 
shocks, and 3) proper time linearization concepts. 

The objective of the present paper is to pre- 
sent a numerical treatment of the unsteady full 
potential equation that properly takes into ac- 
count the importance of the above three listed 
items. The paper discusses a local time linear- 
ization procedure for treating the time derivative 
term, a flux biasing concept based on sonic condi- 
tions (instead of the usual density biasing proce- 
dures) for the treatment of spatial derivative 
terms, and a split boundary condition procedure 
for incorporation into an approximately factored 
implicit algorithm. The resulting unsteady method 
has application not only to unsteady problems in 
transonics, but also to generate the starting 
blunt body solution for a full potential super- 
sonic marching code. 6 The unsteady method of this 
paper, when combined with the steady methods of 
Refs. 4-6, provides a complete treatment of the 
full potential equations. 


Results are presented for flows over cylin- 
ders, spheres and airfoils at various Mach num- 
bers, and some comparisons are made with Euler 
solutions. Use of the split boundary condition 
technique, combined with the flux biasing con- 
cepts, 12,13 has produced a very robust method 
which, even for a difficult case with a fish-tail 
shock at the trailing edge of an airfoil, did not 
require any user specified ‘’constants", such as 
the ones discussed in Ref. 2. 


The paper also presents the results from a 
"hybrid" calculation, wherein the spherical blunt 
body solution from the unsteady code has been 
effectively used to provide the initial data plane 
for a supersonic marching calculating performed 
over the Shuttle Orbiter geometry. 


II. Formulation 


The two-dimensional/ ax i symmetric unsteady full 
potential equation written in a body-fitted coor- 
dinate system represented by t * t, £ ■ C(x,y,t) 
and n * n(x,y,t) takes the form 

( j\ + (p TV + (o 7>„ + 30 Jj “ 0 (l) 

where. 


S * 0 for two dimensions, * 1 for axisymmetric 


p * density * [1- 


U * C + a 


+ n-D) 

n 

+ a. 


1/(Y-1) 


c -ire ' “ 12 V 0 ” 0 + c t 
V \ + a 12*t: + a 22V 9 “ v + T, t 


‘a 


-2 r 2 

C x + V a 12 


c n + r T! ; 

xx y y 


, 2 . 2 
a 22 * * *y 

J * Jacobian * C n - C r\ 
x y y y 


The density 0 and the fluxes oU and pV are 
complicated nonlinear functions of <b , the velocity 
potential. Hence, to solve for from Eq. (1) 
will require a local linearization. 


Let ' n f be the running index in the time di- 
rection, 'k' in the C direction and ’j’ in the n 
direction leading out of the surface. The objec- 
tive is to solve Eq. (1) for at the current 

time plane, knowing the information at n, n-i, 
n-2,... planes. 


• Manager, Computational Fluid Dynamics Group, Associate Fellow AIAA. 

Copyright © American Institute of Aeronautics and 
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A. Treatment 


of (^) in Eq. (1) 


/£nH +1 m 

bx \J ; 


over airfoils or blunt objects), p* is always set 
to p and requires no internal iterations at n+1 
plane. 

The only unknown in Eq. (5) is A4>- 


<a 1 -«b 1 >{<£) n+1 - <§>”} -0b 1 l(5) n -(^) n - 1 l 
a^t^b^ATj + At 2 ) 


where 


ai 

bl 

e 

e 

Ati 
At 2 


■ (Ati + At 2 ) 2 

• Ati* 

* 0 for first order time accuracy 
- 1 for second order accuracy 

n+1 n 
■ t - T 


In order to write Eq. (2) in terms of $ n+1 , a 
local time linearization procedure is introduced. 


C. Treatment of 7— (p 7) in Eq. (1) 

nT} J 


, V.n+1 . $ V n+ \ 

(p J ) n ’ ^ (o * -=T> 




,n+l ; j+l/2,k 


+ a 22 (A* +*”l T1 )^ +1 / 2>k 


( 6 ) 


~ * 


Similarly, P 4 + 1/2 k a modified density based on 
flux biasing. * 


p n+l i p " + (|£) n a* + ... (3) 

where “ (A - * n )» The term (^) is a dif- 
ferential operator given by 


/3fiA n 

K br 


- £_ (2- + 
2 uT 

a 






(4) 


where ’a' is the local speed of sound. 

Combining Eqs. (3) and (4), the nonlinear 
density function in terms of 4> has been linear- 
ized, and the coefficients multiplying A<t> are 
evaluated at the known previous time level. To 
maintain conservation form, both pn+l and p n 
appearing in the first term of Eq. (2) are lin- 
earized according to Eq. (3). 


B. Treatment of 77 (p 7) in Eq. (1) 
oc, J 

TT n +l 

, U.n+l 5_ , p 1 ,k+l/2 u .1,k+l/2 . 

(P ‘ * ( J jf W/2 


ac (a iif A6 + 


11 1 


(5) 


+ a 12 |'A6 + ♦ n } r ,)l J>k+ i /2 

where p n+ * . . appearing in the spatial derivative 
j,k+l/Z 

terra has been linearized to T *' e s y m * , ol 

~ appearing over p denotes that the density has 
been modified to produce the necessary artificial 
viscosity. The modified density is obtained from 
a flux biasing concept to be described later in 
this paper. For a genuine unsteady problem (where 
a time asymptotic steady^state does not exist), 
initially, p* is set to p and then^subsequently 
iterated to convergence by setting p* to the pre- 
viously iterated value of p at the current n+1 
time plane. For problems where the steady state 
exists and is of interest (steady transonic flow 


D. Flux Biasing Procedure 

The density p appearing in Eqs. (5) and (6) is 
modified according to a flux biasing formula given 
by 


for U > 0 


fpq - AC 


j,k+l/2 1 j>k+1 / 2 


(pi) 

AC r j ,k+l/2 


for U < 0 


for V > 0 


q j,k+l/2 


fpq + ACA (pq) 1 j >k+ i/ 2 

oC 


*j+l/2,k “ q j+1/2>k foq " AT1 At, ^J+l/2,k 


for V < 0 






where 


(pq)~ * (pq-p*q*) if q > q* 

* 0 if q < q* 

$ * backward difference operator 
% * forward difference opertor. 


The quantities p*q*, p* and q* represent sonic 
values of the flux, density and total velocity, 
respectively. These sonic conditions are given by 
(using the density and speed of sound relation- 
ships) 


(q*) 


2 1 + ll ^ I M 2 (1-2* t ) 

^M 2 


( 8 ) 


P* 


(q* tO 


2/<v-l) 
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Note that for steady flows, the sonic condi- 
tions p* and q* are only a function of the free- 
stream Mach number, and for a given flow they are 
constants. For unsteady flows, p* and q* need to 
be computed everywhere due to the presence of 4> x 
in Eq. (8). 


Now, the following four cases can be defined. 

a) Subsonic Flow (q < q* at (j,k+l/2) and 

(j, k-1/2)) for U > 0, the modified density in Eq. 

(7) becomes 

^j,k+l/2 ij >k+1 / 2 0<1 l' k+1/2 

[(oq) j,k+l/2 " (pq) j, k-1/2 1 * (9 

- p, k+1 /2 (since (pq)~ at (j, k+1/2) and 
(j,k-l/2) iszero, according to Eq. (8). 


b) Supersonic Flow (q > q* at (j,k+l/2) and 
(j,k-l/2) for U > 0, 


p j, k+1/2 “ 

J. 

q j,k+l/2 


1 ^ pql j, k+1/2 


Upq-p*q *) j(k+1/2 


- (pq~P*q*) j )k -i/2 ^ 


(pq) i.k-l/2 

q j,k+l/2 


+ 


1 

q j, k+1/2 


[ tip *q*)j ( k+l /2 


- (p * q * ) j ,k-l/2 11 


(L0) 


For steady supersonic flows where (p*q*) is a 
constant everywhere, Eq. (10) reduces to 


p j,k+l/2 


0 j » k— 1/2 


| q j ,k 1 /2 | 
q j,k+l/2 


(ID 


c) Transition Through Sonic Line (q > q* at k+1/2 
and q < q* at k-1/2. Refer to Fig. la). For U > 
0, 

q j.k+l/2 , k+1/2 ^ P<1 

- [(pq-p*q*> j>k+l/2 - ( /^,k-i/2 ]1 (12) 

p *q* 

q j ,k+l/2 

d) Transition Through Shock (q > a* at k-1/2 and 
q < q* at k+1/2. Refer to Fig. lb). For 0 > 0, 


/SONIC LINE 

/ 

E / H 
/ 

k-1/2/ k+1/2 k + 3/2 

-• X -j—% X • X •- 

k-1 / k k+1 k + 2 


a) TRANSITION THROUGH SONIC LINE 



b) TRANSITION THROUGH SHOCK 


Fig. 1 Notation for flux biasing. 


p j( k+i /2 o - i j>k+l/2 f(oq) j( k+i/2 

' k+1/2 " <O q ~P * q * ) j,k-l/2 

(13) 

■Pj.k+1/2 + i j)k+1/2 ,D ^* qM j.k-l/2 


For steady flows where p*q* is a constant, it 
can be shown that at a pure supersonic point (case 
b above), the flux biasing procedure and the usual 
density biasing technique of Ref. 2 are identical. 
Using the following relationships: 



(14) 


*£. 

Pic 


is. ia » s_ 

?iq ?)C 2 

a 


q q. 


05) 


fc (pq) 


(p q-p *q*) 

(pq) (for steady flows only) 

OL 


- p r q + pq^ 


(16) 


Using Eq. (15) 


Se (pq) ' 


P c q fl - V 


P r q (1 - 2 

M 


■I (17) 
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Using Eq. (17), thus, for a pure supersonic point, 
Eq. (7) becomes, for U > 0 


and R consists of various known terms at n, n-1 
and n-2 levels. 


P j.k+1/2 


q j,k+l/2 


f p q “ AC P C q (1 " 2 }t j,k+l/2 

(18) 


Using v 


(1 - — r-), Eq. (18) can be written as 
M 


Equation (20) is solved at the (n+1) plane in 
two steps. 


Lg A4> - R, Step 1 
L A4> ■ A4> , Step 2 




+ A 


j,k 


p j,k+l/2 “ p j,k+l/2 " v(p j,k+l/2 
” (1 ‘ v) p j,k+l/2 + vp j,k-l/2 ‘ 


P j,k-l/2 ) 

(19) 


Equation (19) is the usual density biasing tech- 
nique of Holst. 2 However, while transitioning 
through a sonic line (case c) or through a shock 
(case d), the flux biasing procedure of Eq. (7) 
accurately monitors the sonic conditions p* and 
q*, as given by Eqs. (12) and (13). 

The advantages of flux biasing over the 
density biasing* scheme are: 

1. Does not require any user specified constants 
(the parameter 'c* in Ref. 2) that depend on 
the severity of the flow. 

2. Provides a monotone profile through the shock 
wave (for details, see Refs. 12 and 13). 

3. Allows for larger Courant numbers to be taken 
during the calculation (by not producing unde- 
sired pressure overshoots at shocks, which 
could cause instability during transient 
calculations) . 

4. Provides a two point transition through a 
shock wave. 


Both L and L result in tridiagonal matrices. 

C T) 

F. Boundary Condition 

The direction C is along the body surface and 
n leads out of the surface. The boundary condi- 
tion for the Lg operator usually depends on the 
problem. It can be a plane of symmetry condition, 
a periodic condition or a jump in <!> condition that 
goes with a lifting airfoil. All these conditions 
are easily implemented in the Eg operator. Spe- 
cial attention needs to be given for the boundary 
condition that is required in the L^ operator. 

For inviscid flows, the surface flow tangency con- 
dition dictates that the contravariant velocity, 

V, be zero at the body. Implementation of the 

condition, V - 0, in the operator is a crucial 
step in achieving a true implicit scheme. 

Usually, the boundary condition V - 0 is set only 
in the right hand side term R of Eq. (2), and a 
careless or no boundary condition treatment is 
imposed in the left hand side operator. 2 In 
the present method, the condition V * 0 is imposed 
on both sides of the Eq. (20). Let j * J denote 
the body point. Then, 


V “ (a 12*C + a 22*ipJ,k “ ° 


or 


( Vj,k " " ( a 22 *C ^ J ,k 


(21) 

( 22 ) 


A detailed mathematical description of the 
flux biasing procedure can be found in the works 
of Osher 12 and Hafez. 13 

E. Implicit Approximate Factorization 

Combining the various terms of Eq. (1) given 
by Eqs. (2-6) will result in a fully implicit 
model. This is solved using an approximate fac- 
torization implicit algorithm. After some 
rearrangement of the terms, the factored implicit 
scheme becomes 


L„L A$ * R 
C rt 


( 20 ) 


where 


r , . A T , 5 , a ft p* d i 

[i + ATj u + J a 11 grl 


L n “ [1 +AT i v fe + f fej^ a 22 l^ 1 


6 - - ( 


T n+1 f n v2 j,k 

J (a At ^) J 


a * ( 1 - 9 ) +9 [ {a^-b^ (AT j + At^/A^ )/{ a^t^ } ] 


Now, the following adjustments are made for body 
points (using Eqs. (21) and (22)): 

1. Set V - 0 In Eq. (4). 

** 1 O n 

2. Replace a 12 (** + * n ).n ^ (A* + * ) 

ln Eq. (5). a 22 ' 


3 * Replace (p y) by ^ [*- l a 12 (M + * ) £ 

+ a 22 (A* + *)^]j +1 /2 in Eq ’ (&)• This 
assumes (0 j)j_ 1/2 - -<P 7> J+1 /2* 

The approximate factorization at a body point, 
J, is done after the above three boundary condi- 
tion steps are implemented. This leads to the 
following approximate factored scheme for a body 
point, J. 


L_L A4> 
C T) 


(23) 


where 
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+ « JL i*( a . J1 ) 5_j 

6 8C J '•11 a 22 ; e>C J 


Hi [1 + An 8 a 22 aC^+1/2^ 

R • modified form of R in Eq. (20) after 
imposing V a * 0. 

In Eq. (23), the boundary condition is split 
between the two operators, L and L • Even for 
nonorthogonal mesh at the body (a^ 11 * the 
boundary condition is set implicitly. This allows 
for large time steps, or Courant numbers, to be 
taken during the calculation. Imposing Eq. (21) 
directly into the L^ operator in the approximate 
factorization of Eq. (20) would lead to an ex- 
plicit treatment for the part a^j. * n Eq. (21). 

For all calculations to be reported in this 
paper, the above split boundary condition approach 
has proved to be an efficient and robust method, 
even for highly nonorthogonal grid at the body. 

III. Results 


Results are presented for various steady flows 
computed using the unsteady code for time- 
asymptotic steady state solution. The time-step 
Atj appearing in Eq. (2) is computed from a pre- 
scribed Courant number. When impulsively starting 
the calculation from free-stream, the Courant num- 
ber is usually set at 30 ~ 50, and as the itera- 
tion proceeds, the Courant number is rapidly in- 
creased to very large values such as 500 ~ 2000. 

Figures 2 and 3 show the results of a flow 
over a cylinder at M w * 0.4 and 0.45. At ^ • 

0.4, the flow is barely critical, and a comparison 
with an efficient Euler code 14 is excellent. At 
^ * 0.45, a shock is formed on the cylinder sur- 
face. The flux biasing procedure of Eq. (7) pro- 
duces a shock with no overshoots and requires no 
user specified 'constants’ 2 to stabilize the cal- 
culation. Calculations of Figs. 2 and 3 required 
40 time-step iterations. 

Figures 4 to 6 show results for supersonic 
flows over a sphere at different low Mach numbers 
1.08, 1.2 and 1.4. The density distribution for 
all these cases are compared with a bench mark 
type Euler calculation.^- 5 The present full poten- 
tial calculation required 80 time-steps to con- 
verge for all these cases. It is worth noting at 
this point that the Euler calculations of Ref. 15 
for Mach number 1.08 required in excess of 20,000 
interations. 

Figure 7 shows the flow over NACA 0012 airfoil 
at Mach 0.75 and a * 2°. The pressure distribu- 
tion compares well with other full potential meth- 
ods. 2 Since there are no user specified param- 
eters in the dissipation part of the calculation 
(flux biasing), there is only one flow solution 
for the shock in the present method, and it has no 
overshoots. This calculation required 200 itera- 
tions to converge (residual 10 5 ). 


Figures 8 and 9 show the flow over NACA 0012 
airfoil at « 0.98 and a - 0° and 2°, respec- 
tively. The fish-tail shock pattern reported in 
Ref. 2 is nicely reproduced for the zero degree 
case with 250 time-step iterations. At 4° angle 
of attack, the trailing edge shock pattern is 
changed, as seen in Fig. 9. The trailing edge 
shock on the lower surface is weakened, while the 
one on the upper surface grew In strength compared 
to a m 0° case. Even these severe flow cases with 
complex shock shapes required only minimal time- 
sttep iterations and produced no unnecessary 
'wiggles 1 near shocks. 

Figure 10 shows the schematic of a hybrid cal- 
culation where the unsteady code is used to gen- 
erate the blunt body solution for setting up the 
initial data plane for the full potential marching 
code. 16 A hybrid calculation was performed for 
the flow over the Shuttle Orbiter at ^ ■ 1.4, a ■ 
0°. The results of Fig. 6 were used as a starting 
solution for the marching calculation. The nose 
region geometry and the pressure distribution 
along the leeside symmetry of the Orbiter are 
shown in Figs. 11 and 12. 

IV. Conclusions 


An Implicit method for the unsteady full po- 
tential equation is presented. The method employs 
a local time linearization process, a flux biasing 
technique for generation of artificial viscosity, 
and a split boundary condition scheme consistent 
with the approximate factorization algorithm. Re- 
sults for flows over airfoils, cylinders and 
spheres are presented, along with a 'hybrid' cal- 
culation performed for the Shuttle Orbiter, using 
the unsteady and steady codes. Extensions of this 
work into three dimensions is currently in pro- 
gress. Appropriate relaxaton schemes 14 will re- 
place the approximate factorization procedure in 
three dimensions. 
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MACH NUMBER 



Fig. 2 Mach number distribution for cylinder 
flow at M^ * 0.4. 



Fig. 3 Mach number distribution and contours for 
cylinder flow at M m «* 0.45. 
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Fig. k Density distribution and Mach number contours for flow over a sphere 
at - 1.08. 
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Abstract 

The unsteady form of the full potential equation is 
solved in conservation form using implicit methods based 
on approximate factorization and relaxation schemes. A 
local time linearization for density is introduced to enable 
solution to the equation in terms of <f> } the velocity poten- 
tial. A novel flux biasing technique is applied to generate 
proper forms of the artificial viscosity to treat hyperbolic 
regions with shocks and sonic lines present. The wake is 
properly modeled by accounting not only for jumps in <£, 
but also for jumps in higher derivatives of <f > , obtained from 
requirements of density continuity. The far field is mod- 
eled using the Riemann invariants to simulate nonreflecting 
boundary conditions. Results are presented for flows over 
airfoils, cylinders, and spheres. Comparisons are made with 
available Euler and full potential results. 

I. Introduction 

Nonlinear aerodynamic prediction methods based on 
the steady form of the full potential equation are used 
regularly for treating transonic 1,2 and supersonic 3-6 flows 
over realistic wing-body configurations. Numerical algo- 
rithms to compute the unsteady form of the full poten- 
tial equation are still in a developmental stage, and several 
researchers 7- 11 have recently made significant progress in 
this area. There are several issues involved in the con- 
struction of a robust and efficient numerical algorithm for 
the unsteady full potential equation. They are: I) proper 
treatment of boundary conditions in a nonorthogonal grid 
system, 2) correct formulation for the production of ar- 
tificial viscosity to capture sharp shocks, 3) proper time 
linearization concepts, 4) unsteady wake treatment, and 
5) nonreflecting outer boundary conditions. 

The objective of the present paper is to present a nu- 
merical treatment of the unsteady full potential equation 
that properly takes into account the importance of. the 
above listed items. The paper discusses a local time lin- 
earization procedure for treating the time derivative term, 
a flux biasing concept based on sonic conditions (instead 
of the usual density biasing procedures) for the treatment 
of spatial derivative terms, split boundary condition pro- 
cedures consistent with approximate factorization schemes, 
unsteady wake models with proper jumps in <p and higher 
derivatives of <f> taken into account from density relation- 
ships, and nonreflecting unsteady far field procedures based 
on Riemann invariants derived from the characteristic 
theory. 
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Use of unsteady methods has application not only to 
unsteady problems, but also to time asymptotic steady 
state solutions. If the unsteady method is constructed prop- 
erly (robust and efficient), then it can even be made to 
generate steady state solutions faster than methods based 
on the steady form of the equation. Also, in the unsteady 
method, since the time direction is always present, the hy- 
perbolicity of the unsteady full potential equation will allow 
one to obtain solutions to problems across the Mach number 
range (subsonic, transonic, and supersonic), whether steady 
or unsteady. The unsteady method of this paper, when 
combined with the steady marching method of Refs. 4-6, 
provides a complete treatment of the full potential 
equation. 

Results are presented for flows over cylinders, spheres, 
and airfoils at various Mach numbers, and some compar- 
isons are made with Euler solutions. Use of the split bound- 
ary condition technique, combined with the flux biasing 
concepts 12,13 , has produced a very robust method which, 
even for a difficult case with a fishtail shock at the trail- 
ing edge of an airfoil, did not require any user-specified 
“constants”, such as the ones discussed in Ref. 2. 

The paper also presents the results from a “hybrid” 
calculation, wherein the spherical blunt body solution from 
the unsteady code has been effectively used to provide the 
initial data plane for a supersonic marching calculation per- 
formed over the Shuttle Orbiter geometry. 


II. Formulation 

The two-dimensional/axisymmetric unsteady full po- 
tential equation written in a body-fitted coordinate system 
represented by r = f, f = $(x t y t t) and tj = f?(x,y,f) takes 
the form 


(j), + ('7), + ('7), + s V- 0 (,) 


where 

S = 0 for two dimensions, = 1 for axisymmetric 


p — density = 


1 - (Ur 


+ 0 <t><; + V <t>ri - 1 




U = + a + <*12<^ij ; 0 = U + ft 

V = *}t + <* 12 4><; + <*22^; V = V + T} t 

<*11 fx $y j ftc TJx d* $yf}y 

<*22 = + nl 

J = Jacobian = $ x rj y - ; y rj y 


The density p and the fluxes pU and pV are com- 
plicated nonlinear functions of <p t the velocity potential. 
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Hence, to solve for <j> from Eq. (1) will require a local 
linearization. 

Let ‘n’ be the running index in the time direction, ‘A:’ 
in the <; direction, and in the t] direction leading out of 
the surface. The objective is to solve Eq. (1) for at 

the current time plane, knowing the information at n, n- 1, 
n - 2, ■ - • planes. 

A. Treatment of ^ in Eq. (1) 
a_ / £\ n + l _ 

dr \ j) 

(a, - 6b,) { (f) n+ ‘ - (5)"} - 06. {($)" - (?)“" } 

a t A Ti - Obi (Atj + Ar 2 ) 


where p"£+i /2 appearing in the spatial derivative term has 
been linearized to py.fc+ 1 / 2 * The symbol " appearing over 
p denotes that the density has been modified to produce 
the necessary artificial viscosity. The modified density is 
obtained from a flux biasing concept to be described later 
in this paper. For a genuine unsteady problem (where a 
time asymptotic steady state does not exist), initially, p is 
set to p n and then subsequently iterated to convergence by 
setting p to the pre-time plane. For problems where the 
steady state exists and is of interest (steady transonic flow 
over airfoils or blunt objects), p is always set to p n and 
requires no internal iterations at the n + 1 plane. 

The only unknown in Eq. (5) is A <j>. 


where 

a x = (Ar, + At 2 ) 2 
b x ~ A t\ 

0 = 0 for first order time accuracy 
$ — 1 for second order accuracy 
At 1 = r n + l - r n 
At 2 - r n - r n “ l 


(2) C. Treatme n t of ^ {pj ) 

( V\ n+1 . *d /_V tt+, \ 

V p 7]„ = ^V^J )+i /^ 


= £;{j^( a,2{A * +r} < 


(6) 


+ <*22 { + ^ n } f 




In order to write Eq. 2 in terms of ^ n+1 , a local time 
linearization procedure is introduced. 

where A <p = - <p n ). The term is a differential 

operator given by 



where ‘a’ is the local speed of sound. 


Similarly, py+i/ 2 .ic is a modified density based on flux 
biasing. 

D. Biasing Procedure 

The spatial derivative terms given by Eqs. (5) and (6) 
are central differenced expressions about the node point 
(/, k) and are symmetric operators. For shocked flows and 
for treatment of hyperbolic regions, these operators are 
desymmetrized by introducing the biased value of density 
in the upwind direction. This will create the necessary ar- 
tificial viscosity to form shocks and exclude the expansion 
shock. The biased value of density p can be obtained in 
several ways. Some of them are presented here. 


Combining Eqs. (3) and (4), the nonlinear density 
function in terms of <p has been linearized, and the coef- 
ficients multiplying A<p are evaluated at the known previ- 
ous time level. To maintain conservation form, both p n+l 
and p n appearing in the first term of Eq. (2) are linearized 
according to Eq. (3). 


R. Treatment of (py) 


/ vy+ l £ ( bk + i„ ugUjA 

a f { J,. fe+I/2 J 


+ <*i2 {A <f> + <t> n } ri ^ | 


*+1/2 


( 5 ) 


a) Density Biasing 2 (in the ^-direction) 

a (8P\ 

Pk + 1/2 — Pfc+1/2 T I — J , 

\ 0 Wfc+i/ 2 

(7) 

where M is the local Mach number. For U > 0, the - sign 
and backward differencing («— ) is used in Eq. (7), while for 
U < 0, the + sign and (—>) operator is used. 

b) Directional Flux Biasing 


P 




(8) 
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c) Streamwise Flux Biasing 


p= } 


(9) 


where S is the local streamwise direction. Equation (9) can 
be rearranged as 


_ 1 
P ~ q 


pq *{^l< + ^ ri 'k} ( ' >qr } m 


where Q = \/U 1 2 4* V 2 . 

In Eqs. (8) and (10), the term (p?)”" is defined to be 


(pq) = pq - pV 
= 0 


if q>q* 
if q<q* 


( 11 ) 


For steady supersonic flows where (/?*<?*) is a constant 
everywhere, Eq. (14) reduces to 


f?y.fc-i/2l 

Pj\k+if2 - Pj\k- 1/2) (• 

1 <h\k+if2 J 


(15) 


3) Transition Through Sonic Line (q > q* at k - h 1/2 
and q < q* at k — 1/2. Refer to Fig. la.) For £/ > 0, 


^M+i/2 = 


1 


Sj.*+l /2 

- ((/*? - P’q'hk+i/2 ~ \f 4 )lk-ii 2 \) 
p’q 


{{pq)i.k+ 1/2 


(16) 


1 /2 


The quantities p*q% p*, and <7* represent sonic values 
of the flux, density, and total velocity, respectively. These 
sonic conditions are given by (using the density and speed 
of sound relationships) 


(<n 2 


1 4- (1 - 2<ftr - - W0,,) 

2±iM 2 

2 iKi oo 


P' = (g*A/ 0O ) 2/,1 - l> . 


( 12 ) 


4) Transition Through Shock (7 > g* at k - 1/2 and 
q < q" at 4- 1/2. Refer to Fig. lb.) For U > 0, 

P).k+ 1/2 = {(P?)j,Jt+l/2 

9 j.fc+1/2 

- lp^)Jfc+l/2 - (w - PVW 1 / 2 I} (17) 

= Py,fc+i/2 d {p<? ~ P*?*}y,fc~i/2 

<lj\k+l /2 


Note that for steady flows, the sonic conditions p * and 
q * are only a function of the freestream Mach number, and 
for a given flow they are constants. For unsteady flows, p* 
and q* need to be computed everywhere due to the presence 
of <p T and other unsteady terms in Eq. (12). 

The density biasing based on flux, Eq. (10), is more 
accurate than the one presented in Eq. (7), since it is based 
on sonic reference conditions. To illustrate the flux biasing 
procedure for various situations, Eq. (8) is considered. 


1) Subsonic Flow (q < q* at (j, k + 1/2) and (;, k - 

1/2)) for U > 0, the modified density in Eq. (8) becomes 


hk+i/2 = «y.*+i/a{ (M)i ’* +l/a 

- [(pq) j.k +l /2 («Wl/ a ]} 


(13) 


= P].k+ 1/2 (since (pq)~ at {j,k+ 1/2) and (;, k - 1/2) is 
zero, according to Eq. (11), 

2) Supersor ic Flow {q > q * at (;, k + 1/2) and ( j,k - 
1/2)) for U > 0, 


Pj.k + l /2 = 


-{(P<l)j,k+ 1/2 


<?/,*+ 1/2 

~[(pq — p‘q‘)j.k+i/2 — {pq — p* )y.jt— 1/2I } 
(pq)j.k- 1/2 


(14) 


1/2 1/2 

-(p’q'U- 1 /2 }]. 


K(pYWi/2 


For steady flows where p‘q * is a constant, it can be 
shown that at a pure supersonic point (case 2 above), the 


/ SONIC LINE 

/ 

E / H 
/ 

k-1/2/ k+1/2 k + 3/2 

-• x f • x • x • 

k-1 / k k+1 k + 2 


a) TRANSITION THROUGH SONIC LINE 



b) TRANSITION THROUGH SHOCK 


Fig. 1. Notation for flux biasing. 
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flux biasing procedure, Eq. (8), and the usual density bias- 
ing technique of Ref. 2, Eq. (7), are identical. 



dp _ dp dq _ p 
d$ dq d$ a 

Q 

= — (p<j) (for steady flows only) 
= Pf? + M< 


(18) 

(19) 


( 20 ) 


Using Eq. (20) 



( 21 ) 


Using Eq. (21), thus, for a pure supersonic point, Eq. (8) 
becomes, for U > 0 


1 

Pj.k+ 1/2 — 

£>,*+ 1/2 

Using v - (l - ), Eq. (22) can be written as 

Py.A+i/2 = Pj\k+ii2 “ ^(p/Jfe+i/a — Py,fc-i/a) 

= (1 - V)pj,k+lf 2 + ^ Pj,k— 1 / 2 * 


( 22 ) 


(23) 


Equation (23) is the usua] density biasing technique of 
Holst 1 2 3 4 . However, while transitioning through a sonic line 
(case 3) or through a shock (case 4), the flux biasing pro- 
cedure of Eq. (8) accurately monitors the sonic conditions 
p* and q * , as given by Eqs. (16) and (17). 

The advantages of flux biasing over the density 
biasing 2 scheme are: 

1. Does not require any user-specified constants (the pa- 
rameter ‘c’ in Ref. 2) that depend on the severity of 
the flow. 

2. Provides a monotone profile through the shock wave 
(for details, see Refs. 12 and 13). 

3. Allows for larger Courant numbers to be taken during 
the calculation (by not producing undesired pressure 
overshoots at shocks, which could cause instability dur- 
ing transient calculations). 

4. Provides a two point transition through a shock wave. 
A detailed mathematical description of the flux bias- 
ing procedure can be found in the works of Osher 12 and 
Hafez 13 . 


E. Unsteady Wake Treatment 

Figure 2 shows the schematic of a wake cut behind the 
trailing edge of an airfoil. This wake cut has to be properly 
modeled in the unsteady formulation. An expression for 
the jump in <p across the wake cut can be derived by re- 
quiring that the pressure be continuous across the cut. In 
the full potential framework, this results in the continuity 
of density. Equating the density p u = pi at any point along 
the wake (refer to Fig. 2), one can write 

2r t 4* (U u + Ut) {fa) u - U t T< = 0 (24) 

assuming ^ 0, where [/] stands for the jump in the 

quantity given by (/„ - f t ). Equation (24) is integrated 
from the trailing edge to the downstream boundary (along 
TE in Fig. 2) to obtain the T distribution along the wake 
cut. To maintain stability, the £ derivatives in Eq. (24) 
are upwind differenced. For a steady flow, Eq. (24) will 
result in a constant value for T along the wake given by 
r ss <j>j> — 4> t , in Fig. 2. 

Beside the T evaluation, solution to the unsteady equa- 
tion also requires information on (<£ n ) at a wake point. Re- 
ferring to Fig. 2 for notation, one can write the following 
using Taylor’s series expansion. 

<t>2 = fa - (fa)u + (^ijn )«/2 H 

fa ~ 4>z ~ H — 

The subscripts ‘u’ and *£’ stand for upper and lower, re- 
spectively. 

For the chosen coordinate system, requiring that 
(<Mu = and defining fa - fa = T, using Eqs. (25) 

one can write 

W,). = + |^|/2} (2C) 


Equation (26) requires an estimate for the jump in <p nri at 
the wake cut. This can be obtained by setting the jump in 
the equation to be zero at a wake point. 


( 5 ). ♦ (n), + ( 4 ) 



(27) 



Fig. 2. Wake and outer boundary treatment. 
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Writing down Eq. (27) in terms of the transformation met- 
rics and requiring that p u = pi , the following relationship 
is obtained. 


\<Pm\ = - 


P & 23 


(28) 


For steady flows, the jump in <f> nrJ along the wake is zero 
because F f — 0. For unsteady flows, to maintain accuracy, 
Eqs. (24)— (23) have to be employed. 


F. Far Field Boundary Condition 

Along the outer boundary ABODE in Fig. 2, appro- 
priate Riemann invariants are prescribed. The concept can 
be explained by considering the Cartesian form of the full 
potential equation cast in terms of the Riemann invariants. 



Equation (29) implies that along the (u + a) positive char- 
acteristics the quantity is invariant, and along 

the (u-a) negative characteristics the quantity 
is invariant, known as the Riemann invariants. Usually, 
along the outer boundary, the Riemann invariant that cor- 
responds to the positive characteristics with respect to the 
inward normal can be prescribed as a boundary condition. 

For an arbitrary coordinate system (r, such as the 
one in Fig. 2, the following boundary conditions are 
appropriate. 

U 2 

__ _| - a = constant along AB 

^/a^ 7-1 

U 2 

7 = + -a — constant along ED (30) 

v^rr 7-1 v ’ 

V 2 

7 = -i -a = constant along BCD 

V ^22 7 -1 


Equation (30) is nonlinear in nature. Hence, to implement 
the Riemann invariant boundary condition, a linearization 
technique similar to Eq. (3) is employed. For example, 
equating the right hand side of Eq. (30) to the freestream 
value, along BCD one can write 


( a ‘ 2 ^ 


2 £7 + a 22 




-T-{h +v §;+ v ii) ** 




(31) 

( V 2 \" 

freestream V + 7 ~ 1 7 


Riemann invariant boundary conditions is better than pre- 
scribing 0oo from the compressible vortex solution, and will 
serve as a nonreflecting boundary condition. 

G. Relaxation and Approximate Factorization Schemes 

When all the terms of Eq. (1) are put together, it will 
be in terms of the unknown A <j>. It can be written as 

F(M)+R{r,t n - l , -) = o , M = r +l -r. (32) 

The Newton iterative procedure for solving Eq. (32) be- 
comes 

dF 

— (At n+l - A^) = -*(**, *»- 1 , ■ ■ •) - F( A*‘) (33) 


The off-diagonal elements of which cannot be accom- 
modated within the tridiagonal setup of SLOR, can be set 
to zero. In Eq. (33), A<£* represents the intermediate it- 
erative value of A <fi. For steady state problems, can 
be set to A<£ n . If all the off-diagonal terms of are set 
to zero, then the relaxation process becomes the point Ja- 
cobi iteration. Results based on the relaxation procedure 
of Eq. (32) are presented in this paper. 

Another procedure to solve Eq. (32) is the approximate 
factorization technique. This can be written as 

L <; L rt A<f> = R (34) 


L,= 


Lr, - 


where 

3 a 9 p' 9 

+ - — ya u — j 

3 a dp' 9 

i + Ari K_ + ___a 22 _j 

P = ~ (j»+‘(a»Ar l ) a )y. fc 

a = (1 - 9) + 5[{o, - 6 i(At! + Ar a )/Ari} / (a t - 6 t }] 
Equation (34) is solved at the (n + 1) plane in two steps. 


A0 = R , Step 1 
Lr, A<j> = A<^ , Step 2 

*% l = tik + ^ 

Both and L v result in tridiagonal matrices. 


H. Body Boundary Condition 

For inviscid flows, the surface flow tangency condition 
dictates that the contravariant velocity, V ) be zero at the 
body. Implementation of the condition, V = 0, in the 
Lr, operator is a crucial step in achieving a true implicit 
scheme. Usually, the boundary condition V = 0 is set only 
in the right hand side term R } and a careless or no bound- 
ary condition treatment is imposed in the left hand side L n 
operator 2 . In the present method, the condition V = 0 is 
imposed on both sides of Eq. (34). Let j = J denote the 
body point. Then, 

V = ~ ® (25) 


The finite differenced form of Eq. (31) will provide an es- 
timate for (A0) along the outer boundary. Use of the 
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or 


(36) 



Using Eqs. (35) and (36), and the relationship 



(37) 


the and the L n operators in Eq. (34) can be modified to 
the form, for a body point J, 


L f — R 


(38) 


where 


L< 


x + ^ u Ts 


a'd 

0 


Ln = 


'LL (a 

J V 11 “22 ) 

2 a ( p d \ 

&V 0 \J a 22 ds) 


J+ 1/2J 


In Eq. (38), the boundary condition is split between the two 
operators, and L n . Even for ncnorthogonal mesh at the 
body (a i2 ^ 0), the boundary condition is set implicitly. 
This allows for large time steps, or Courant numbers, to be 
taken during the calculation. 




Fig. 3. Mach number distribution for cylinder flow at M = 
0.4. 


III. Results 

Computer codes based on both the relaxation method 
(Eq. (33)) and the approximate factorization method 
(Eq. (34)) have been constructed for two-dimensional and 
axisymmetric problems. The grid around the geometry can 
be either a C-type (Fig. 2) or an O-type. At present, the 
relaxation method is somewhat slower (at least 50%) in 
convergence than the approximate factorization code. How- 
ever, the future implementation of multigrid techniques 17 
and implicit relaxation concepts 14 to the present relaxation 
code can make it competitive to approximate factorization 
methods in three-dimensional applications, where the ap- 
proximate factorization methods with triple factorization 
can be less flexible to handle complex geometries 14 . 

The unsteady code has been applied, at present, only 
to steady state problems. Calculations involving unsteady 
motions such as plunge, pitch, and oscillating flaps are cur- 
rently in progress. For steady state problems, the time step 
A Tj appearing in Eq. (34) is computed from a prescribed 
Courant number, usually set much greater than one (~ 50). 

Figure 3 shows the result for a flow over a cylinder at 
Moo =0.4. The flow is barely critical, and a comparison 
with an efficient Euler code 14 is excellent. Figures 4 and 5 
show results for supersonic flows over a sphere at low Mach 
numbers of 1.08 and 1.4. The density distribution for these 
cases are compared with benchmark Euler calculations 15 . 
The present full potential code required approximately 80 
time steps to converge (residual < 10" 5 ). It is worth noting 
at this point that the Euler code of Ref. 15 requires in 
excess of 20,000 iterations to perform the low Mach number 
calculation of 1.08. 

Figure 6 shows the pressure distribution obtained over 
the NACA 0012 airfoil at Moo = 0.8 and a = 0° with 


a grid of 84 points around the airfoil and 18 in the re- 
direction. The comparison with an Euler code 14 is good. 
Figure 7 shows the results for M Q 0 = 0.75 and a = 2° 
over the same airfoil. Even with a crude grid, the for- 
mation of a strong shock without any overshoots is made 
possible by the use of flux biasing concepts. Calculations 
of this type require no user-specified “constants” to in- 
crease or decrease the amount of dissipation. Depending 
on the strength of the shock, the flux biasing automati- 
cally chooses the right amount of dissipation, since it is 
based on sonic reference conditions. The perfect matching 
of pressure contours across the wake cut (Fig. 7) illustrates 
the correctness of the unsteady wake model described in 
this paper. Figure 8 shows some difficult cases with fishtail 
shocks. 

The unsteady code can also be effectively used to gen- 
erate the blunt body solution, the outflow of which is to be 
prescribed as an initial data plane for a full potential su- 
personic marching code 6 . Figure 9 shows the schematic of 
such a hybrid calculation. The flow over the entire Shuttle 
Orbiter with a blunt nose has been simulated at Mo© = 1.4 
and a = 0°. The results of Fig. 5 were used as a starting 
solution for the marching calculation 16 . The nose region 
geometry and the pressure distribution along the leeside 
symmetry of the Orbiter are shown in Figs. 10 and 11. 

Simulation of unsteady phenomena, such as flutter and 
control surface oscillations will be presented in the future. 

TV. Conclusions 

A computational treatment for the unsteady full po- 
tential equation is presented. The method employs a local 
time linearization, flux biasing concepts for generation of 
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artificial viscosity, unsteady wake treatment, outer bound- 
ary conditions based on the Riemann invariants, and re- 
laxation and approximate factorization algorithms. Use of 
the code for problems with steady state solution has been 
very effective and computationally fast. Extensions cf this 
work to simulate unsteady phenomena such a3 flutter, and 
to three dimensions to treat wings, and wing-body combi- 
nations, are currently in progress. 
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